---
title: "Simulation study to assess the properties of the t-interval"
author: ''
output:
  html_document:
    toc: yes
  pdf_document: null
  pdf_presentation: null
  slidy_presentation:
    css: slidy.css
    widescreen: yes
  word_document: null
---


## Preliminary explanations

Let us start today's activity with a few words of explanation.

In statistics, we often want to represent data by a <b>Model</b>

Doing this allows us, to a certain extent, to say something that goes beyond the data itself.

To carry out this exercise, we select a model in our bag of models.


Most models have parameters, which characterize the shape of their probability distribution.

Hence, if I select the normal distribution as my model (as we will do here), I need 2 parameters, the mean and the variance of the normal distribution, to represent my data.

But we do not know these parameters.

Therefore we estimate these parameters, assuming that the model is correct (that the data was sampled from a normal distribution).


Having estimated the mean and the variance, we can calculate 'the error bar' of our estimator of the mean.

The standard method for doing this is called t-confidence interval for the mean estimator of a normal distribution.

So, in the exercise we propose today, we will use R to make simulations and check, empirically, if the analytical theory is correct.


We will also play a trick. We will generate data from another model, and see how wrong our confidence interval is.


To check the quality of our confidence interval we are going to estimate its 'coverage'.

The idea is simple: repeat several estimations of the mean, using synthetic=simulated data and see if the confidence intervals does indeed contain the true (this time known!) value of the parameter. 

We will collect the proportion of replicates in which the calculated interval indeed contains the mean. This is the empirical coverage of the t-confidence interval.


And now, on with calculations in R.







## Carrying out a simulation study

  The generic framework for carrying out a simulation study is

```{r,eval=F}
  set #Nbatch#, the number of simulation batches Nbatch
  set #n#, the sample size of each simulation batch
    for (i in 1: Nbatch){ 
       draw a sample of size n from a distribution
       do something using the sample
       }
   aggregate results across simulation batches
```


## Example: simulation study to evaluate the coverage of the t-interval

The t-interval is given by

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

where $t_{\alpha/2,n-1}$ cuts off an area $\alpha/2$ in the upper tail of the
$t$ distribution with $n-1$ degrees of freedom.

### when samples are drawn from a normal population

When the data are a sample from a normal population with mean $\mu$ and
standard deviation $\sigma$, then $100(1-\alpha)$% of t-intervals constructed will
contain the true but unknown population mean $\mu$.  The **nominal**, or named, level of
the t-interval is $100(1-\alpha)$%.

The following simulation study draws 100 samples of size 10 from a normal
population with mean 3, standard deviation .5.  For each sample, a 95\%
confidence interval for the mean is calculated.  The 100 intervals are
plotted, and a count is made as to how many of the intervals contain
the true mean.  The proportion of intervals containing the true mean
is the "empirical coverage" of the procedure.

```{r}
 plot(0,0,xlim=c(2,4),ylim=c(0,100),type="n",
       ylab="simulation batch", xlab="when sampling from normal population")
 Nbatch=100 #fix the number of simulation batches
 n=10       #fix the sample size for each batch
 truemean=3; sd=.5; alpha=.05 #fix some parameters for this specific simulation
 cover=0

for (i in 1:Nbatch){
   data=rnorm(n,truemean,sd) #draw a sample
   #now do something with the sample              
   tint=mean(data)+c(-1,1)*qt(1-alpha/2,n-1)*sd(data)/sqrt(n)
   count1=ifelse(tint[1]<truemean&&truemean<tint[2],1,0)
   lines(tint,c(i,i))
   cover=cover+count1
 }

lines(c(truemean,truemean),c(0,100))
 empirical=cover/Nbatch 
 print(paste("empirical coverage is ",100*empirical,"%"))
```

The estimated, or **empirical**, coverage of the t-interval is
`r empirical`%. If we did more simulation batches the empirical coverage would
typically be closer to the nominal coverage.  This is a property of the t-interval
when sampling from a normal population.

### when samples are drawn from a non-normal population

What if we sample from a non-normal popluation?
The following does the similar simulation, but where samples are
drawn from  $\chi^2_2$, the $\chi^2$ distribution with 2 degrees of freedom.  This is a very
non-normal population, as can be seen from the following.

```{r}
x=seq(0,5,length.out=200)
plot(x,dchisq(x,2),type="l",ylab="f(x)",main="chi-squared_2 density function")
```

Now carry out the same simulation study, but with samples drawn from the $\chi^2_2$
distribution.

```{r}
 plot(0,0,xlim=c(-1,6),ylim=c(0,100),type="n",
     ylab="simulation batch", xlab="sampling from chi-squared distribution")
 truemean=2; n=10; Nbatch=100; alpha=.05; df=2
 cover=0
 for (i in 1:Nbatch){
   data=rchisq(n,df=2)
   tint=mean(data)+c(-1,1)*qt(1-alpha/2,n-1)*sd(data)/sqrt(n)
   lines(tint,c(i,i))
   cover=cover+ifelse(tint[1]<truemean&&truemean<tint[2],1,0)
 }
  lines(c(truemean,truemean),c(0,100))
 empirical=cover/Nbatch  
 print(paste("empirical coverage is ",100*empirical,"%"))
```
 

With $n$ moderately large, the empirical coverage will be fairly close to
the nominal coverage.
This is because of the central limit theorem, which says that the distribution
of the sample mean is quite close to a normal distribution even if the
underlying population is non-normal.  For smallish values of $n$, such as
was used here, the empirical and nominal coverage can be quite different.


# Hints

## Dissecting the code

You may be completely lost with the first cell of R code

Therefore, I have dissected it for you here:

```{r}

# first initialize the plot
 plot(0,0,xlim=c(2,4),ylim=c(0,100),type="n",
      ylab="simulation batch", xlab="when sampling from normal population")


# we will call 

 Nbatch=100 #fix the number of simulation batches
 n=10       #fix the sample size for each batch
 truemean=3; sd=.5; alpha=.05 #fix some parameters for this specific simulation
 cover=0

for (i in 1:Nbatch){
  
# draw n sample from the normal distribution
  
   data=rnorm(n,truemean,truesd) #draw a sample

# estimate the parameters
   estmean = mean(data)
   estsd   =   sd(data)
   
# calculate the t-confidence interval for the estimated mean:   
   tint(1)=estmean + tn*estsd/sqrt(n)
   tint(2)=estmean - tn*estsd/sqrt(n)
   
# now count1 is 1 if this segment contains the true parameter (true mean) or else zero
   
   count1=ifelse(tint[1]<truemean&&truemean<tint[2],1,0)
   
# increment the counter : cover counts hits
   cover=cover+count1
   
# add graphics: draw the segment a bit below
   lines(tint,c(i,i)) 
   
 }

# draw a vertical line where the true mean is
lines(c(truemean,truemean),c(0,100))
# all plots are now finished


# calculate the coverage probability
 empirical=cover/Nbatch 
# spit out the result:
 print(paste("empirical coverage is ",100*empirical,"%"))

# punchline?
# is the estimated coverage the same as the announced coverage?
 

 
```
 
 

## Wrapping the code into a function


```{r}
bete <- function(n=n,alpha=alpha,Nbatch=Nbatch,
                 truemean=truemean,truesd=truesd,showplot=showplot){

# first initialize the plot
if(showplot == 1){
 plot(0,0,xlim=c(2,4),ylim=c(0,Nbatch),type="n",
      ylab="simulation batch", xlab="when sampling from normal population")
}

 #fix some parameters for this specific simulation
 cover=0

for (i in 1:Nbatch){
  
# draw n sample from the normal distribution
  
   data=rnorm(n,truemean,truesd) #draw a sample

# estimate the parameters
   estmean = mean(data)
   estsd   =   sd(data)
   
# calculate the t-confidence interval for the estimated mean: 
   
   tnx <- qt(1-alpha/2,n-1)
   tint = estmean + c(-1,1)*tnx*estsd/sqrt(n)
   
# now count1 is 1 if this segment contains the true parameter (true mean) or else zero
   
   count1=ifelse(tint[1]<truemean&&truemean<tint[2],1,0)
   
# increment the counter : cover counts hits
   cover=cover+count1
   
# add graphics: draw the segment a bit below
   if(showplot == 1){
   lines(tint,c(i,i)) 
   }
   
 }

# draw a vertical line where the true mean is
 if(showplot == 1){
lines(c(truemean,truemean),c(0,100))
 }
# all plots are now finished


# calculate the coverage probability
 empirical=cover/Nbatch 
# spit out the result:
 print(paste("empirical coverage is ",100*empirical,"%"))

# punchline?
# is the estimated coverage the same as the announced coverage?
 
 theoretical = 1. - alpha
 print(paste("theoretical coverage is ",100*theoretical,"%"))
 return (empirical) 
}

# we will call 

 Nbatch=100 #fix the number of simulation batches
 n=10       #fix the sample size for each batch
 truemean=3; 
 truesd=.5; 
 alpha=.05;
 showplot=0;
 
bete(n=n,alpha=alpha,Nbatch=Nbatch,truemean=truemean,truesd=truesd,showplot=showplot) 
 
```
 
 OK
 
 Now can you write a function that calls bete 1000 times?
 
 What do you conclude?
 
 
