---
title: "R Lecture "
output: html_notebook
---


Discrete distributions: 
Let's plot the probability mass function of some discrete distributions
```{r}
x<-NULL
for (i in 1:20){
  x[i]<-dbinom(i,20,0.5)
}
plot(x)
```


```{r}
x<-NULL
for (i in 1:20){
  x[i]<-dgeom(i,0.5)
}
plot(x)
```


```{r}
x<-NULL
for (i in 1:50){
  x[i]<-dpois(i,5)
}
plot(x)
```



```{r}
x<-NULL
for (i in 1:50){
  x[i]<-dnbinom(i,20,0.5)
}
plot(x)
```



We can also plot the cdf's by prefixing with a p: 

```{r}
x<-NULL
for (i in 1:20){
  x[i]<-pbinom(i,20,0.5)
}
plot(x)
```

Exercise: Plot the cdf's of some other discrete random variables. 


Continuous Distributions: 

The common continuous distributions are norm, exp, gamma, unif, beta, t, chisq, f, cauchy. 

The prefix "p" gives the value of the cdf at a point. For example: 

```{r}
pnorm(0,0,1) #the z-score of 0 is 50%
```

```{r}
pt(3, 8) #The probability that t<=3 where t is a t-distributed random variable with 8 degrees of freedom. 
```


Let's plot the pdf of  the T distributions. We can do this by sampling the distribution with the R prefix, and then plotting a histogram: 



```{r}
x<-rt(10000,8) #create a vector with 1000 random samples of a t-distribution with 8 degrees of freedom.
hist(x,breaks=100, prob=TRUE) #plot a histogram of these values
lines(density(x),lwd=3, col="red") #Adds a density function based on the data. 
```


```{r}
x<-rexp(10000,0.01) #create a vector with 1000 random samples of an exponential random variable with lambda = 0.01.
hist(x,breaks=100,prob=TRUE) #plot a histogram of these values
lines(density(x),lwd=3, col="red")
```


```{r}
x<-rchisq(20000,25) #create a vector with 2000 random samples of an chi^2 random variable with 25 degrees of freedom.
hist(x,breaks=150, prob=TRUE) #plot a histogram of these values
lines(density(x),lwd=3, col="red")
```





Quantile Function: 

The quantile function of a distribution is the inverse of the CDF. We have already used this impicitly when we find critical values. For example, given some alpha, z_alpha is the quantile function of the standard normal evaluated at 1-alpha. 

The quantile function uses the prefix "q". 

For example: 

```{r}
qnorm(.95) #give the critical value z such that P(Z<z)=0.95
```

```{r}
qt(0.9, 7) #give the critical value t_{0.1,7} such that P(t_7<t_{0.1,7})=0.9
```


Simulation Studies

Simulation studies have been used in past lectures to illustrate various statistical phenomenon such as the CLT. 

Simulation studies require us to sample a distribution randomly many times over. 

We use the "r" prefix to do this. The "r" prefix returns an IID sample of a random variable with a given size. 

For example: 

```{r}
rnorm(10,mean=5,sd=3)# generates a random sample of size 10 from a normal distribution with mean 5 and sd = 3
```

Sometimes we want a lot of samples. There are several ways to do this. One way is with a for loop: 

```{r}
means <- rep(NA, 100) #Creates a blank vector of length 100. 
for (i in 1:1000) {
  x <- rnorm(10, 5, 3) #for loop creates 1000 instances of an iid sample of                            size 10
  means[i] <- mean(x)  #stores the mean of each sample. 
}
```

```{r}
hist(means,breaks = 50)
```

Another way is to just create a matrix, where we consider each row to be a sample: 

```{r}
# Generate 10000 samples of size 10 from the normal population
Samples <- matrix(rnorm(10*10000, 5, 3), ncol=10) #ncol specifies the number of columns, so we have 10 columns and 1000 rows, i.e. 1000 samples of size 10. 
```

```{r}
xbar <- rowMeans(Samples) #takes the mean of each row and stores them in a vector xbar. 
```



Let's use this to do a simulation study of CI's 

```{r}
#decide on a
a<-0.05 #a=0.05 is a 95% CI
# Calculate the 1000 a% confidence intervals
CI.a <- cbind(xbar-qnorm(1-(a/2))*3/sqrt(10), xbar+qnorm(1-(a/2))*3/sqrt(10))
#remark: the cbind function is just applying the argument to each element of xbar (which is a vector) and returning a vector of pairs (the lower and upper endpoints of the interval)
```
Notice that "qnorm(1-a/2)" is the quantile function that gives the value of z_{\alpha/2} for alpha=a i.e. a two sided 95% CI on the mean of a normal distribution.


```{r}
qnorm(1-(a/2))
```


```{r}
# Look at the first 10 intervals
CI.a[1:10,]
```


Now let's count how many times the true mean of our population (mu=5) actually is in our CI: 

```{r}
# Count how many include the true population mean mu=5
sum(CI.a[,1]<5 & CI.a[,2]>5)
```