---
title: "R Programming Basics - Useful Builtin Functions for Statistics"
output:
  html_document:
    toc: yes
  pdf_document: null
  pdf_presentation: null
  slidy_presentation:
    css: slidy.css
    widescreen: yes
  word_document: null
---


# Vectorized Arithmetic - most arithmetic operations in R work on vectors.

Here are a few commonly used summary statistics. 

Please try to use them, e.g. to each column of a dataframe.


```{r}
testVect = c(1,3,5,2,9,10,7,8,6)

min(testVect) # minimum
max(testVect) # maximum
mean(testVect) # mean
median(testVect) # median
var(testVect) #variance
sd(testVect) # standard deviation
# HEY YOU CAN GENERATE VECTORS BY selecting the column of datasets !
# if cars is a dataframe with a column named speed :
vect1 = cars$speed
# another example
vect2 = cars$dist

# now you can treat each vector as an empirical distribution
# and apply your statistical functions:
quantile(vect1) # sample quantiles 
summary(vect1)  # same thing, plus the mean

cov(vect1,vect2) # sample covariance between vect1 and vect2
cor(vect1,vect2) # sample correlation coefficient
```

# Calculations on probability distributions - area under curve, percentiles, random sampling

## distributions
    - *binom* : Binomial Distribution (Pascal)
    - *pois*  : Poisson Distribution (Poisson)
    - *unif* :  Uniform Distribution 
    - *exp* : Exponential Distribution 
    - *norm* : Normal Distribution (Gauss, Laplace, Legendre)
    - *chisq* : Chi-Squared Distribution (Pearson)
    - *t*   : t Distribution (Student)
    - *f*  : F Distribution (Fisher)
    
    
## Why do we mean by distributions?

Load the course summary (main index) in your browser and go to
the 'Supplementary material section'. 

Load and read <b>Hey, let's talk about distributions!</b>

This material is for fun. Nothing there is required for the course.

    
## key prefix letters
    # *d* : probability density function/probability mass function. 
    # *p* : cumulative distribution function. 
    # *q* : quantile function. Returns percentiles, also known as quantiles.  
    # *r* : random number generation


## Explaining pdf, cdf and so on

<ul>

<li><b>q (quantile):</b> What is the value of the OUTCOME (observed quantity) 
             that 'corresponds' to a GIVEN probability (input for quantile)
             
for example: what is the 80% quantile of canadian men heights?

If $X$ is the random value of the observation (i.e. the measured height of a randomly sampled canadian man), we are looking for (a R function that gives) the value of the height $H$ such that $P(X < H) = 0.80$. What is your guess for $H$? How would you estimate it? 


EXAMPLE

Assume that $X \sim N(1.7, 0.5)$

How do I use R to calculate the value $q$ such that $P(X < q)=0.8$?

```{r}
sdval=sqrt(0.5)
qnorm(0.8,mean=1.7,sd=sdval)
```

Is this a credible result for human height? Why?


<li><b>p (probability):</b> This is the reverse of quantiles. We give the outcome (i.e. the height) and compute the PROBABILITY THAT THE OUTCOME is less then H.
For example, we give the height: $H=1.781m$ and we want an algorithm (R function) that calculates $P(X < 1.781cm)$. In 2014, this was about $0.5$. The probability that a random variable is less than the given numerical value $x$ is called the <emph>cumulative distribution function</emph>: $F(x)=P(X < x)$.


EXAMPLE:

Say $X \sim Norm(0,1)$

How do I use R to calculate $P(X < 1.5)$?

```{r}
pnorm(1.5)
```

Geometric interpretation (AUC):

<img src="http://daad.wb.tu-harburg.de/uploads/pics/cdf-pdf-1_03.gif">


<li><b>d (probability density function):</b> This is is the derivative of the cumulative density function. This function $f(x)$ can be thought of as the probability of realizing an outcome in the immediate vicinity of the numerical value $x$, i.e. $f(x)dx = P(X \in [x,x+dx])$.

EXAMPLE:

Suppose that $X \sim Binom(N=4,p=0.3)$

What is the probability that $X=3$?


```{r}
dbinom(3,size=4,prob=0.3)
```

How would you check that this is correct?`
             
<li><b>r (random sample):</b> Here we are looking for a way to PRODUCE A NEW OBSERVATION, i.e. a new height for another (randomly sampled) man. We must choose some model (a generative distribution model) that represents the relative importance of men heights in the canadian male population, and must find an algorithm to sample (i.e. produce a new example of output) from this distribution. R does this for us, for a couple of well-known distributions.

EXAMPLE:
<a href="https://chemicalstatistician.wordpress.com/2013/04/14/checking-the-goodness-of-fit-of-the-poisson-distribution-for-alpha-decay-by-americium-241/"></a>

In my house, I have a smoke detector.

This detector uses the radioactive decay of a small amount of Americium (Nova-Scotium to be more exact) to detect the presence of smoke (which disrupts the electrical contact caused by radioactive decay in a small ion chamber)

The decay reaction (just for the curious) is:

241Am   →   237Np   +   4He


Assuming that alpha particles are produced at an <b>average</b> rate of 0.8392 emissions by second, and that the count of emissions in intervals of 5 seconds obey a Poisson distribution, how would I ask R to produce a sequence of random counts that mimics the production of alpha particles in my smoke detector for 15 intervals?

Can you try to use the rpoisson function to do this?

```{r}
#rpois(...)
```

This exercise is called : random sampling from a given distribution (here: a Poisson distribution).


</ul>

## HEY BUT CAN R FUNCTIONS REALLY TELL ME THE PROBABILITY OF EVENTS in my dataset?

No.

The R functions make ideal calculations based on mathematical models.

This may differ from the 'truth' because in many cases, we do not know from what distributions the data came from.

However, probability theory and statistics provide ways to estimate these distributions and the targets are often the model distributions mentioned here.

Let us take an example.

Suppose that $X$ is the height of a canadian male in 2014. Will R tell me that $P(X < 1.781cm)=0.5$? Absolutely not. R does not know the distribution. You will have to tell R what kind of distribution to use. If you pretend (very good approximation, see live histograms or the work of Sir Francis Galton on human stature) that $X$ comes from or follows a normal distribution, then <b>you will have to provide the value of the parameters: mean and variance of the distribution</b>. Finally, you will just ask R to calculate $dnorm$ for those parameters.

## HOW CLOSE are real distributions from model distributions?

For an answer, see the supplementary material (Let's talk about distributions!)

## Examples of manipulation of the statistical distributions in R

   - What is the area under the standard normal curve to the left of 1.5?

```{r}
pnorm(1.5)
```

   - What point cuts off an area .975 to its left under the standard normal curve?  this point also cuts off an area .025 to its right under the standard normal curve.

```{r}
qnorm(.975)
```

   - What is the area under the standard normal curve to the left of 1.96?  To the left of -1.96?
  What is the area under the curve betwen -1.96 and 1.96?  

```{r}
pnorm(1.96)
pnorm(-1.96)
pnorm(1.96)-pnorm(-1.96)
```

  
   - Simulate 100 observations from the normal distribution with mean $\mu=3$ and standard deviation
$\sigma=1.5$.  Calculate some summary statistics for the data, and also show a normal QQ plot,
which can be used to assess normality.


```{r}
data=rnorm(100,mean=3,sd=1.5)
summary(data)
qqnorm(data)
qqline(data)
```

   - Plot the standard normal density function $\phi(z)$ from $z=-3.1$ to $z=3.1$, where

$$\phi(z)=\frac{1}{\sqrt{2 \pi}} e^{-.5 z^2}$$

```{r}
z=seq(from=-3.5, to=3.5, length.out=300)
dz=dnorm(z)
plot(z,dz,xlab="z",ylab="phi(z)",main="standard normal curve",type="l")
```

   - What is the probability that a normal random variable with mean 4 and standard deviation 1.5 is greater than 2.5?

```{r}
1-pnorm(2.5,4,1.5)
```

   - What is the probability that a $t$ random variable with 12 degrees of freedom is less than 2?
Greater than 2?

```{r}
pt(2,12)
1-pt(2,12)
```

   - What are the upper and lower 2.5'th percentile of a $t$ random variable with $6$ degrees of freedom?

```{r}
qt(c(.025,.975),6)
```

   - Simulate 200 observations from the uniform distribution on the interval (0,5).  Make a 
histogram of the data with 10 bars.

```{r}
data=runif(200,0,5)
hist(data,nclass=10)
```

   - What about with 1000 observations?  Does it look more uniform?  Make a probability histogram, which
means the bars are normalized so that the total area is 1.

```{r}
data=runif(1000,0,5)
hist(data,nclass=10,prob=T)
```

   - What about with 10000?  The law of large numbers indicates that the 
histogram will look close to uniform for a sufficiently large sample size.

```{r}
data=runif(10000,0,5)
hist(data,nclass=10,prob=T)
```

   - What is the probability that a binomial random variable with parameters $n$ and $p$ is equal to $x$?

$$p(x) = \left ( \overset{n}{x} \right ) p^x (1-p)^{n-x}$$

eg. with $n=10$, $p=.2$, $x=2$

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

   - What is the probability that a binomial random variable $X$ with parameters $n$ and $p$ is less than or equal to $x$?

$$P(X \le x ) = \sum_{k=0}^x \left ( \overset{n}{x} \right ) p^x (1-p)^{n-x}$$

eg. with $n=10$, $p=.2$, $x=2$

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

   - Simulate 20 observations from the binomial distribution with $n=10$ and $p=.2$.

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

## When generating random numbers, use  *set.seed()* to reproduce results.

```{r}
set.seed(100)
rnorm(5)
rnorm(5)

set.seed(100) # reproduce the results
rnorm(5)
rnorm(5)
```


#  Useful functions for working on matrices.


## apply - apply a function over rows or columns of a matrix

   - syntax: apply(matrixname,dimension,FUN)

   - example: for the airquality data set, count the number
   of missing values by row (for each record), and by column (for each variable)
   

```{r}
nmiss=function(x){sum(is.na(x))}
apply(airquality,1,nmiss)
apply(airquality,2,nmiss)

```{r}
mydata=read.csv("http://www-bcf.usc.edu/~gareth/ISL/Auto.csv")

## %*% -inner product.  dot product. matrix multiplication.

```{r}
M1=matrix(1:9,byrow=T,ncol=3)
M2=matrix(1:15,byrow=T,ncol=5)
M1
M2
M1%*%M2
```

## outer - calculate an outer product between two arrays

    outer(X, Y, FUN = "*", ...)
     
    X, Y: First and second arguments for function ‘FUN’.  Typically a
          vector or array.

FUN: a function to use on the outer products.


The function FUN will be called using each element of X and each element
of Y.

Example:  The following calculates a table of values x^y, for each x in
1:10, and each value of y in -2:2.  Row and column names are added to
the table.

```{r}
x=1:10
y=-2:2
powertable=outer(x,y,FUN="^")
rownames(powertable)=x
colnames(powertable)=y
```






<!-- ## some packages (for fun) -->

<!-- - *sqldf*  -->
<!-- ```{r} -->
<!-- library(sqldf) -->
<!-- library(ISLR) -->

<!-- sqlCars = sqldf(" select name,origin,year,cylinders,horsepower -->
<!--                   from Auto  -->
<!--                   where cylinders = 8 and horsepower > 200 -->
<!--                   order by horsepower -->
<!--                 ") -->

<!-- sqlCars -->
<!-- ``` -->

<!-- - *quantmod* -->
<!-- ```{r} -->
<!-- library(quantmod) -->
<!-- getSymbols("AAPL",src="yahoo",from = Sys.Date()-500) -->
<!-- tail(AAPL) -->
<!-- barChart(AAPL) -->
<!-- getQuote("AAPL",src="yahoo") -->
<!-- ``` -->


