###################################################
# Computer Lab 4 – F71SM
###################################################
## In the tasks below, fill in any missing code as
## required.
## Make sure you use the help menus in R (e.g.”?mean”
## will open up a help window for “mean()”).
## In some cases, the numerical answers you will find
## may be slightly different from those in the tutorial
## answers, due to rounding in intermediate steps.
###################################################
# Sampling distn for X.bar – Normal population
# See also notes, Ch6 p.3
# Generate 500 samples of size n=30 from N(10,4) population
# For each sample compute the sample mean
# Plot a histogram of the sample means
# Compute the mean and variance of the sampling distribution
# and compare them with their theoretical values
nsim = 500 # No. of simulated samples
n = 30 # Sample size
xbar.30 = numeric(nsim) # Initialise vector
for (i in 1:nsim){
x = rnorm(n,10,2)
xbar.30[i] = mean(x)
hist(xbar.30,xlim=c(5,15),freq=F)
# freq=F used to produce bars with proportional areas
# this will be useful later
mean(xbar.30); var(xbar.30)
# Repeat code above (incuding for loop) to see the randomness
# in the sampling distribution
# Repeat the exercise in part (i) for sample size n=100
# Add the theoretical pdf of the distribution of Xbar on the histogram
# Remember, this should be a N(10, 4/100) distribution
hist(xbar.100,xlim=c(5,15),freq=F,main=””) # Draw histogram again
title(main=”Histogram of Xbar and pdf of N(10,0.04)”)
x.grid = seq(9,11,length=200) # Grid of values on ax axis
lines(x.grid,dnorm(x.grid,10,sqrt(0.04)),col=”red”,lwd=2) # Draw pdf
###################################################
# Sampling distn for X.bar – Chi-squared population
# See also notes, Ch6 p.6
# Repeat part (ii) from Task 1, by generating 500 samples
# of size n=100 from Chi.Squared(4) population
# Repeat part (iii) from Task 1
# What is the distribution of Xbar now?
###################################################
# Mean Square Error
# See also examples 7.7, 7.17
# Generate 10^4 samples of size n=6 from a Uniform (0,theta) distribution, with theta=2.
m = matrix(0,nrow=10^4,ncol=n) # Initialise matrix, each sample in a row
for (i in 1:10^4){
m[i,] = runif(n,0,theta) # Generate samples in rows
m[1:10,] # Check first 10 samples
# We have seen in lectures (7.7, 7.17) that two estimators for theta are:
# (a) 2*Xbar and (b) X(n) (the largest observation)
# Calculate the 2 estimates from each of the 10^4 samples
theta.a = 2*apply(m,1,mean) # Estimates (a); Find out about ?apply
theta.a[1:10] # Check estimates of theta in first 10 samples
theta.b = … # Estimates (b)
theta.b[1:10]
# Calculate the mean and bias of each estimator and comment.
# Hint: the expected value of the estimator can be approximated
# by the mean of the estimates from the 10^4 samples above
# Calculate the variance and MSE of each estimator and comment.
# Hint: again, the variance of the estimator can be approximated
# by the variance of the estimates from the 10^4 samples above
###################################################
# Method of moments estimation
# See also example 7.9
# The following sample is considered to come from a Gamma(a,b) distribution:
y = c(3.9,3.4,5.2,2.6,3.4,2.7,3.1,2.4,4.2,1.7,2.4,3.3,2.3,2.0,3.4,2.9,2.5,3.9,2.3,1.1)
# Calculate estimates of the parameters a and b using the method of moments.
###################################################
# Maximum likelihood estimation
# See also example 7.6
# Consider X ~ Binomial(10,theta) and observed value x=7
# Write two functions in R calculating the log-likelihood function and
# the likelihood function for given theta, n (no.of trials) and x
# Calculate the log-likelihood and likelihood when theta=0.5, n=10, x=7
# Check your answers visually with graphs in Ch7, p6
ll.f <- function(theta,n,x){ k1 = log(choose(n,x)) log.likel = k1 + x*log(theta) + (n-x)*log(1-theta) return(log.likel) l.f <- function ... # Plot the likelihood and log-likelihhod functions # over a suitable range of theta when n = 10, x = 7 # Check your answers visually with graphs in Ch7, p6 par(mfrow=c(1,2)) # Plotting window layout: 1x2 theta.grid = seq(0.05,1,length=200) # Grid of theta values on x axis # Obtain the MLE by maximising numerically the log-likelihood # Use the optim() function in R # Compare the results to the analytical answer in 7.6 init.value = 0.5 # Need to provide an initial "guess" for theta optim(init.value,ll.f,n=10,x=7,control=list(fnscale=-1))$par # ?optim for help # Ignore the warning message # Also, compare the max loglikelihood value with the graph # and the value calculated with ll.f() at the MLE optim(init.value,ll.f,n=10,x=7,control=list(fnscale=-1))$value