##### Demonstration Script #2a
##### MATH322
#####
##### An investigation into coverage
##### 


### Preamble

# Here is a function that returns TRUE/FALSE,
# depending on whether the given value (x) is
# between the other two (lo and up).
#
# This starts to illustrate how to write our
# own functions in R.
#

isBetween <- function(x,lo,up) { 
  ib = x>lo & x<up
return(ib)
}





### Examples for the concept of coverage:
#   If these examples are slow, then you may want
#   to change the number of iterations. Fewer will
#   make things run faster. 1e4 is really the fewest
#   you will want to use. More gives a more-precise
#   estimate of the coverage.




### Example 1. For mu; X~Normal, sigma known
alpha = 0.05
n = 10

# Distribution parameters
mu = 0
sigma=1

# Initialize variables
covered = numeric()
lo = numeric()
up = numeric()

# Begin Monte Carlo simulation
for(i in 1:1e4) {
 x = rnorm(n, m=mu, s=sigma)
 lo[i] = mean(x) + qnorm(alpha/2)*sigma/sqrt(n)  ## lcl
 up[i] = mean(x) - qnorm(alpha/2)*sigma/sqrt(n)  ## ucl
covered[i] = isBetween(0,lo[i],up[i])
}

# Coverage estimate
mean(covered)
1-alpha





### What is happening?
#   Let us look at this visually... at least for the first several

plot.new()
plot.window( xlim=c(-2,2), ylim=c(0,200))
axis(1, at=0, label=expression(mu), cex.lab=1.5)
abline(v=mu)

for( i in 1:210) {
	segments(lo[i],i, up[i],i, col=2-covered[i])
Sys.sleep(0.25)
}

# Actual test that our coverage is what we would expect
binom.test(x=sum(covered),n=1e4)
# If this confidence interval includes 1-alpha, then
# there is no evidence that these calculations are wrong










### Example 2a. For mu; X~Normal, sigma unknown

alpha = 0.05
n = 10

mu = 0
sigma=1

covered = numeric()
lo = numeric()
up = numeric()

for(i in 1:1e4) {
 x = rnorm(n, m=mu, s=sigma)
 lo[i] = mean(x) + qnorm(alpha/2)*sd(x)/sqrt(n)    #### Oops! This is wrong (but easy)
 up[i] = mean(x) - qnorm(alpha/2)*sd(x)/sqrt(n)    #### What could possibly go wrong???
covered[i] = isBetween(0,lo[i],up[i])
}

mean(covered)
1-alpha

binom.test(x=sum(covered),n=1e4)
# If this confidence interval includes 1-alpha, then
# there is no evidence that these calculations are wrong.
#
# This confidence interval does NOT contain 1-alpha. THUS, 
# this method to calculate confidence intervals for mu
# under these circumstances is wrong!!!
#
# What happened??? Yeppers. It was the Guinness Effect!
#









### Example 2b. For mu; X~Normal, sigma unknown

alpha = 0.05
n = 10

mu = 0
sigma=1

covered = numeric()
lo = numeric()
up = numeric()

for(i in 1:1e4) {
 x = rnorm(n, m=mu, s=sigma)
 lo[i] = mean(x) + qt(alpha/2,  df=n-1)*sd(x)/sqrt(n)    #### Fixed!
 up[i] = mean(x) + qt(1-alpha/2,df=n-1)*sd(x)/sqrt(n)    #### Fixed! =)
covered[i] = isBetween(0,lo[i],up[i])
}

mean(covered)
1-alpha

binom.test(x=sum(covered),n=1e4)

# And so, the interval based on the z-distribution (Ex2) is not 
# appropriate if you have to use the variance of the data in lieu
# of that of the population. William Sealy Gosset discovered this
# back in the early Century XX (dos equis??). 
#
# Here is a citation for his article. We will be looking at it 
# later in this course (if all goes according to plan):
# "The probable error of a mean." (1908) Biometrika 6(1): 1--25.
#







### Example 3. For mu; X~Exponential

# To see what an Exponential distribution
# looks like, here is a graphic of its pdf

a=2

x = seq(0,10,length=1e4)
y = dexp(x, rate=a)
plot(x,y,type="l")
Ex = 1/a
# It is highly non-Normal. So, let us see if the 
# t-procedure works in this extreme case.

alpha = 0.05
n = 10
mu = 0
sigma=1

covered = numeric()
lo = numeric()
up = numeric()

for(i in 1:1e4) {
 x = rexp(n, rate=a)
 lo[i] = mean(x) + qt(alpha/2,df=n-1)*sd(x)/sqrt(n)
 up[i] = mean(x) - qt(alpha/2,df=n-1)*sd(x)/sqrt(n)
covered[i] = isBetween(Ex,lo[i],up[i])
}

mean(covered)
1-alpha

binom.test(x=sum(covered),n=1e4)
# Does the confidence interval reported here contain 1-alpha?
# Nope.
# Is it close?
# It depends on what you are doing with your estimates.
#
















### Example 4. For sigma^2; X~Normal

alpha = 0.05
n = 10

mu = 0
sigma=1

Ex = mu
Vx = sigma^2

covered = numeric()
lo = numeric()
up = numeric()

for(i in 1:1e4) {
 x = rnorm(n, m=mu,s=sigma)
 lo[i] = (n-1)*var(x)/qchisq(1-alpha/2,df=n-1)
 up[i] = (n-1)*var(x)/qchisq(alpha/2,df=n-1)
covered[i] = isBetween(Vx,lo[i],up[i])
}

mean(covered)
1-alpha

binom.test(x=sum(covered),n=1e4)













### Example 4. For sigma^2; X~Beta(a,b)
a=2
b=4

# This is what the distribution looks like
x=seq(0,1,length=1e4)
y=dbeta(x, shape1=a,shape2=b)
plot(x,y,type="l")
# The Beta distribution is used to model proportions.
# It is very flexible. The standard Uniform is 
# equivalent to Beta(1,1).


# The mean and variance of Beta(a,b)
Ex = a/(a+b)
Vx = a*b/( (a+b)^2 * (a+b+1) )

alpha = 0.05
n = 10

covered = numeric()
lo = numeric()
up = numeric()

for(i in 1:1e4) {
 x = rbeta(n, shape1=a,shape2=b)
 lo[i] = (n-1)*var(x)/qchisq(1-alpha/2,df=n-1)
 up[i] = (n-1)*var(x)/qchisq(alpha/2,df=n-1)
covered[i] = isBetween(Vx,lo[i],up[i])
}

mean(covered)
1-alpha

binom.test(x=sum(covered),n=1e4)

# Oops, the coverage is too high in this example. In fact,
# the coverage will ("almost certainly") be too high, regardless
# of a and b, simply because the Beta is bounded between 0 
# and 1.
#
# Why should boundedness matter?













### Example 5. For theta; X~Unif(0,theta)
#
# We know what this distribution looks like, so there
# is no need to see it

theta = 29

x=seq(0,theta,length=1e4)
y=dunif(x, min=0,max=theta)
plot(x,y,type="l")

##
alpha = 0.05
n = 10

covered = numeric()
lo = numeric()
up = numeric()

for(i in 1:1e4) {
 x = runif(n, min=0, max=theta)
 U = max(x)
 lo[i] = (1-alpha/2)^(-1/n)*U
 up[i] = (alpha/2)^(-1/n)*U
covered[i] = isBetween(theta,lo[i],up[i])
}

mean(covered)
1-alpha

binom.test(x=sum(covered),n=1e4)

# How did we do? The confidence interval reported by this last
# line should include 1-alpha. Why? Because we used mathematics
# to get this result.
#
# BUT, you objet, we had another estimator of theta, it was the MOM
# estimator: 2xbar. Is there a confidence interval based on that
# estimator???????
#
# Exactly? No. (Not at this level. The calculations are complicated.)
# Approximately? Yes.
#
# What is that approximation based on? 
# The CLT!!!  =)
#




### Example 6. For theta; X~Unif(0, theta)

alpha = 0.05
n = 10

theta = 29


covered = numeric()
lo = numeric()
up = numeric()

for(i in 1:1e4) {
 x = runif(n, min=0,max=theta)
 lo[i] = 2*( mean(x) + qt(alpha/2,  df=n-1)*sd(x)/sqrt(n) )
 up[i] = 2*( mean(x) + qt(1-alpha/2,df=n-1)*sd(x)/sqrt(n) )
covered[i] = isBetween(theta,lo[i],up[i])
}

mean(covered)
1-alpha

binom.test(x=sum(covered),n=1e4)

# And here, we see the power of the t-test! 
# Even though the data are not generated from
# a Normal process, we were able to estimate 
# the population mean very well.
#
# Since we could estimate the population mean, 
# we were able to estimate confidence intervals
# for it. AND, since we could estimate confidence
# intervals for the population mean, we could 
# estimate linear functions of it. Most specifically,
# we could estimate 2mu, which is also theta, the
# population parameter of interest.
#