---
title: "Practice session: "
author: ''
output:
  html_document:
    toc: yes
  pdf_document: null
  pdf_presentation: null
  slidy_presentation:
    css: slidy.css
    widescreen: yes
  word_document: null
---


## Outline

- Vectors, matrices, tables
- Flow Control
- Built-in statistical functions
- Write functions

#vectors matrices tables


what is v?

```{r}
a=3
b=9
v  = a+c(-1,2,4)*b
```


what does this do?

```{r}
x=c(T,T,F,F)
y=c(T,F,T,F)
M=cbind(x,y,x&y,x|y,!x)
dimnames(M)[[2]]=c("x","y","x&y","x|y","!x")
M
```





# Flow Control

## Conditional execution

```{r eval=T}
fulltime=NA
nclasses=1
if (nclasses >= 3) {fulltime=TRUE
} else {fulltime=FALSE}
fulltime
```


### ifelse

#### *ifelse* is a compact version of *if ... else ...*
#### Syntax
```{r eval = F}
#ifelse(condition, statement1, statement2)
```
There is a subtle difference in that *ifelse* works on vectors.

#### Examples

```{r}
#find the component by component minimum of vectors a and b
a=rnorm(10); 
b=rnorm(10); 
vecminab=ifelse(a<b,a,b); 
print(cbind(a,b,vecminab))
```


if number if negative write 'negative'

if number is positive write '>0'


## Loops


```{r}
W=matrix(1:16,byrow=T,ncol=4)
print(W)

fmat=function(W){
n=nrow(W)


# these loops will generate 16=4x4 states for the matrix
#
    
      for (i in 1:n){
      for (j in 1:n){
        
        W[j,i]=W[i,j]+W[j,i]
        
      #  pass i=1
        #11  = 11+11 # 11 will never be seen again    
        #21  = 12+21 # 21 will be seen again
        #31  = 13+31
        #41  = 14+41
        
      #  pass i=2
        #12  = 21+12 # will never be seen again 
        #22  = 22+22
        #32  = 23+32
        #42  = 24+42 
        
      #  pass i=3
        #13  = 31+13 # 
        #23  = 32+23
        #33  = 33+33
        #43  = 34+43
        
      #  pass i=4
        #14  = 41+14 # 
        #24  = 42+24
        #34  = 43+34
        #44  = 44+44 
        
      }
      }
    
      return(W)
}

print(fmat(W))
```



```{r}



```


# Statistical functions


<ul>
<li>vectorized arithmetics
<li>pdf,auc,percentiles,random sampling
<li>functions of matrices
<li>binom,pois,exp,unif...
<li>d,p,q,r
</ul>


<ul>
<li>empirical coverage
<li>simulation
<li>t-interval
<li>plot
</ul>


what does this code do?


```{r}
print(sample(c(2,5,3), size=3, replace=FALSE))
print(sample(c(2,5,3), size=3, replace=TRUE))

```


what does this code do?


```{r}
pbinom(2, size=10, prob=.2)
```

what does this code do?


```{r}
data=rbinom(2, size=10, prob=.2)
data
```


<a href="https://homepage.cs.uiowa.edu/~jblang/probability.calculators/t.table.htm">t-table</a>
<img src="https://homepage.cs.uiowa.edu/~jblang/probability.calculators/t.tabl2.gif">

$t_{\alpha,k}$

quantiles

```{r}
qt(.975,23)
```



what does this code do?


```{r}
set.seed(87612345)
mydata=rnorm(20,mean=4,sd=.75)
t.test(mydata,conf.level=.95)
```

<ul>

<li>what is the assumption made for data?


<li>are these assumptions justified?


<li>what is a confidence interval?


<li>I used the expression: 'is the rat in the box'. What did I mean?


<li>what is the coverage of a confidence interval?


<li>what does empirical coverage of a CI mean?


<li>what is alpha? what is dof?


<li>what is the formula that R is using in the function t.test?


<li>where is a quantile used? quantile of what? what does a quantile mean?


<li>how could I have a smaller confidence interval?


<li>what would happen if the data does not come from a normal distribution?


<li>Can you give me examples?


<li>but how do I know the distribution of the data?


</ul>


what does this code do?


```{r}
set.seed(117426)
data=rchisq(42,7)
```



<b>what is the assumption made for data?</b>


data is i.i.d normally distributed.

Each data point is a random number 

$$x_i \sim N(\mu=4,\sigma^2 = (.75)^2)$$

I will use the notation 

$$x_i \sim f$$ 
to say the data is i.i.d. fro, distribution $f$

Please note that in general, distribution have 'knobs' (controls or parameters), so a better notation could be:

$$x_i \sim f(.,\theta)$$ 

where $theta$ is a vector of parameters for the distrubution.

For example, a multivariate normal distribution in 2 dimensions requires 5 scalar parameters.


<b>what is a confidence interval?</b>

A confidence interval for a parameter of a distribution (for example the mean of a normal) is an interval in which we believe that the parameter will be with some given probability.

For example, a $90\%$ confidence interval for the  mean $\mu$ of a random sample (data) is an interval $[a,b}$ such that


$$ \mathbb{P}( \mu \in [a,b])  = 0.9$$



<b>I used the expression: 'is the rat in the box'. What did I mean?</b>


The rat is the unknown parameter. The box is the confidence interval. 


<b>what is the coverage of a confidence interval?</b>


It is the real probability that the interval contains the parameter.
This may differ from the nominal or announced probability.


<b>what does empirical coverage of a CI mean?</b>


It is an <b>estimation</b> of the real coverage probability of the CI.
We get it by simulating data and counting how many times the rat is in the box.



<b>what is alpha? what is dof?</b>


for a $90\%$-CI, alpha is $10\%$. The number of degrees of freedoms depends on the rat. It will not be the same if you try to catch a 'mean' or 'a variance' say. For a mean, the dof is the number of data points minus 1.</b>


<b>what is the formula that R is using in the function t.test?</b>

The formula of Student

When $x_1, x_2, \ldots , x_n$ is a sample from a normal distribution
with unknown mean $\mu$ and unknown variance $\sigma^2$, the level
$100(1-\alpha)$% confidence interval for $\mu$ is given by

$$\bar x \pm t_{1-\alpha/2, n-1} \frac{s}{\sqrt n}$$

where $\bar x$ and $s$ are the sample mean and sample standard deviation
of the data, and $t_{1-\alpha/2, n-1}$ cuts off an area $1-\alpha/2$ to
its left under the $t$ curve with $n-1$ degrees of freedom.



<b>where is a quantile used? quantile of what? what does a quantile mean?</b>

see below

<b>how could I have a smaller confidence interval?</b>

Get more data points


<b>what would happen if the data does not come from a normal distribution?</b>

The coverage of your t-CI will not be the nominal coverage


<b>Can you give me examples?</b>

Draw data from any other distribution : poisson, chisquare, binomial, uniform etc.


<b>but how do I know the distribution of the data?</b>

You try fitting several distributions and keep the distribution with the best fit. BUT...


# Writing functions


STEPS for writing functions:

step 1: formulate the problem

step 2: solve the problem

step 3: implement your algorithm in R

NOTA-BENE: You can build your function on top of powerful functions available in base R or in
libraries!


## compute a series

write a function  that calculates 

$$ f(a,n) = \sum_{i=1}^n \frac{1}{i^a}$$


for $a=2$, this is:

$$(1/1)^2 + (1/2)^2 + (1/3)^2 + (1/4)^2 + \ldots$$

compute $ \sqrt{6 * f(a=2,n=20)}$ 

does this remind you of any well-known number?


My solution

```{r}
sa=function(a,n){
  return(sum(seq(1,n)^(-2)))
}

n=20
print(sqrt( 6.*sa(2,n)))
```



```{r}
n=100000
print(sqrt( 6.*sa(2,n)))
pi
```


some amazing speed effects


```{r}
#for (n in (20:100000)){
#s = sqrt( 6.*sa(2,n))
#if(n%%10000 == 0){print(c(s,s-pi))}
#}
```



```{r}
for (n in (20:100000)){
if(n%%10000 == 0){print(sqrt( 6.*sa(2,n)))}
}
```



## compute the agm

You all know that given two numbers, $x,y$, we can define several types of averages:

$$A(x,y)=1/2(x+y)$$
$$G(x,y)=\sqrt{x y}= (xy)^{1/2}$$
Note that ...

$$H(x,y)=\frac{1}{\frac{1}{x} + \frac{1}{y}}$$


Suppose I start with two real numbers

$$a=a_0,b=b_0$$

Now I define a sequence of points

$$a_{n+1}=A(a_n,b_n), b_{n+1}=G(a_n,b_n)$$

Can you write a function that calculates $a_N,b_N$ for any given rank $N$?


```{r}
#agmf=function(arguments here){
   #
   #return()
#}
```

The arithmetic-geometric mean is defined as the common limit of these two sequences

<a href="https://arxiv.org/pdf/1907.04110.pdf"></a>

FUN: Here is a code found online that plays with the AGM function.

Can you guess what it does?

```{r}

library(numbers)

##  Gauss constant
1 / agm(1, sqrt(2))  # 0.834626841674073

##  Calculate the (elliptic) integral 2/pi \int_0^1 dt / sqrt(1 - t^4)
f <- function(t) 1 / sqrt(1-t^4)
2 / pi * integrate(f, 0, 1)$value
1 / agm(1, sqrt(2))

##  Calculate pi with quadratic convergence (modified AGM)
#   See algorithm 2.1 in Borwein and Borwein
y <- sqrt(sqrt(2))
x <- (y+1/y)/2
p <- 2+sqrt(2)
for (i in 1:6){
  cat(format(p, digits=16), "\n")
  p <- p * (1+x) / (1+y)
  s <- sqrt(x)
  y <- (y*s + 1/s) / (1+y)
  x <- (s+1/s)/2
  }

## Not run: 
##  Calculate pi with arbitrary precision using the Rmpfr package
require("Rmpfr")
vpa <- function(., d = 32) mpfr(., precBits = 4*d)
# Function to compute \pi to d decimal digits accuracy, based on the 
# algebraic-geometric mean, correct digits are doubled in each step.
agm_pi <- function(d) {
    a <- vpa(1, d)
    b <- 1/sqrt(vpa(2, d))
    s <- 1/vpa(4, d)
    p <- 1
    n <- ceiling(log2(d));
    for (k in 1:n) {
        c <- (a+b)/2
        b <- sqrt(a*b)
        s <- s - p * (c-a)^2
        p <- 2 * p
        a <- c
    }
    return(a^2/s)
}
d <- 64
pia <- agm_pi(d)
print(pia, digits = d)
# 3.141592653589793238462643383279502884197169399375105820974944592
# 3.1415926535897932384626433832795028841971693993751058209749445923 exact

## End(Not run)
```


## solve the cubic quation


Write a function that prints one root f a cubic polynomial


<a href="https://math.vanderbilt.edu/schectex/courses/cubic/">link</a>




$$x^{3}+px+q=0$$



$$t=x+{\frac {b}{3a}}$$
$$p={\frac {3ac-b^{2}}{3a^{2}}}$$
$$q= {\frac {2b^{3}-9abc+27a^{2}d}{27a^{3}}}$$
$$4p^{3}+27q^{2}>0$$


$${\sqrt[{3}]{-{\frac {q}{2}}+{\sqrt {{\frac {q^{2}}{4}}+{\frac {p^{3}}{27}}}}}}+{\sqrt[{3}]{-{\frac{q}{2}}-{\sqrt {{\frac {q^{2}}{4}}+{\frac {p^{3}}{27}}}}}}$$

$$ax^3+bx^2+cx+d=0$$

The discriminant of the general cubic 

$$18\,abcd-4\,b^{3}d+b^{2}c^{2}-4\,ac^{3}-27\,a^{2}d^{2}$$


It is the product of 
$$a^{4}$$ 
and the discriminant of the corresponding depressed cubic.


$$\Delta _{0} = b^{2}-3ac $$
$$\Delta _{1} =2b^{3}-9abc+27a^{2}d$$


$$C={\sqrt[{3}]{\frac {\Delta _{1}\pm {\sqrt {\Delta _{1}^{2}-4\Delta _{0}^{3}}}}{2}}}$$



$$x=-{\frac {1}{3a}} (b+C+{\frac {\Delta _{0}}{C}} )$$

Can you write a function that accepts a depressed cubic as input and returns three things:

- a message

- one real root of a cubic polynomial?

- the discriminant



```{r}
#cubicroot=function(arguments here){
   #
   #return()
#}

#discriminant=
#return(NULL)}

```



## approximate an integral


The integral of a function $f(x)$ can be approximated by a Riemann sum

$$\int_{a}^b f(x)dx \approx \sum_{i=1}^n f(x_i) dx_i$$

How would you write a function that calculates the riemann sm at any gvn anl f a fnctn?


$$\int_a^b g'(x)dx = g(b)-g(a)$$

if $f(x)=g'(x)$, $g$ is called a primitive of $f$

What if $g=sin$, what is $f$?

Hence 

$$\int_a^b  sin(x)dx = \cos(b)-\cos(a)$$



```{r}
dumm=function(x,f){
 output=f(x)+2
 return(output)
  }
dumm(pi/2,sin)
```


Implement a function that calculates the riemann sum at rank n of a function of your choice

What output might you want to return?


```{r}
rsum=function(f,n){
# output=f(x)+2
# return(output)
  }
# how would you test your function?
```


```{r}
bread=function(x){
# implement any function of your choice
# return()
  }
# how would you test your function?

```

```{r}
butter=function(x){
if(x == 0 ){return(1)} else{
return(sin(x)/x)
}
}
# how would you test your function?
butter(0.01)

```

what is the integral of butter?


```{r}
rsum=function(f,n){
# output=f(x)+2
# return(output)
  }
# how would you test your function?
```



How about computing the area under an empirical (not analytical) curve?



```{r}

x <- seq(-1., 1., by = 0.05)

#set.seed(100)
y = exp(-x^2/2) / sqrt(2*pi)

library(pracma)
area <- trapz(x, y)
area

```




```{r}

z <- seq(-10., 10., by = 0.05)

xs <- sample(z,20,replace=F)

x <- sort(xs)

#set.seed(100)
y = exp(-x^2/2) / sqrt(2*pi)

library(pracma)
area <- trapz(x, y)
area

```

# Vectorized arithmetics




```{r}
```

