Stat 6501 Homework 3
Enrique Rodriguez
November 4, 2015
6c, 7c, 16c, 35,47abc,50, 52c, 55, 58
6.
Suppose that \(X\sim bin(n,p)\)
- If \(n=10\) and \(X=5\), plot the log likelihood function.
Solution:
Calculating the log-likelihood for the Binomial PMF yields:
\[l(p) = \log \binom{n}{X} + X\log(p) + (n-X)\log(1-p)\] substituting in for the values given we have: \[l(p) = \log(252) + 5\log(p) + 5\log(1-p)\]
7.
Suppose that \(X\) follows the geometric distribution.
- Find the asymptotic variance of the MLE.
Solution
By results in a proven theorem we have that the assymptotic variance is given by:
\[-\frac{1}{El''(p)}\]
So the problem reduces to finding the value of the denominator.
First we can formulate the log-likelihood function to be: \[l(p) = n\ln(p) + (\sum_{i = 1}^nX_i - n)\ln(p)\] Taking the first and then second derivative of this with respect to \(p\) yields: \[l'(p) = \frac{n}{p} - \frac{\sum X_i - n}{1-p} \Rightarrow\] \[l''(p) = -\frac{n}{p^2} - \frac{\sum X_i - n}{(1-p)^2}\] Now taking the expectation of this (we will take the negative and inverse of this at the end): \[E\Big[-\frac{n}{p^2} - \frac{\sum X_i - n}{(1-p)^2}\Big] = -\frac{n}{p^2} - \frac{\sum E[X] - n}{(1-p)^2}\] \[ = - \frac{n}{p^2} - \frac{n(1/p) - n}{(1-p)^2}\] \[n\Big[- \frac{1}{p^2} - \frac{(1/p) - 1}{(1-p)^2}\Big]\] \[= n -\Big[\frac{(1-p)^2 + p - p^2}{(1-p)^2p^2}\Big] =-\frac{n}{(1-p)p^2}\]
Now applying the negative and inversing we get the assymptotic variance to be: \[ \frac{(1-p)p^2}{n}\]
16.
Consider an i.i.d sample of random variables with density function \[f(x|\sigma) = \frac{1}{2\sigma}\exp\Big(-\frac{|x|}{\sigma}\Big)\]
- Find the asymptotic variance of the mle.
Solution
Using what we formulated in the finding the mle we know that the log-likelihood function of the density function is given by:
\[l(\sigma) = -n\ln(2\sigma) - \frac{\sum |X|}{\sigma}\]
Now we find the first and second derivative of this with respect to \(\sigma\)
\[l'(\sigma) = -\frac{n}{\sigma} + \frac{\sum |X|}{\sigma^2}\] \[l''(\sigma) = \frac{n}{\sigma^2} - 2\frac{\sum |X|}{\sigma^3}\]
Now taking the expectation of this (we take the inverse and negative at the end):
\[E\Big[l''(\sigma)\Big] = E\Big[\frac{n}{\sigma^2} - 2\frac{\sum |X|}{\sigma^3}\Big]\] \[ = \frac{n}{\sigma^2} - \frac{2}{\sigma^3}\sum E[|X|] \;\; **\]
We now tackle the problem of finding the exptected value in the equation.
\[E\Big[|X|\Big] = \int_{-\infty}^\infty x\frac{1}{2\sigma}\exp(-\frac{x}{\sigma}) = \frac{1}{2\sigma}\int_{-\infty}^\infty x\exp(-\frac{x}{\sigma})\]
But now taking a close look we can rearrange this so that it looks like something we are familiar with, namely the exponential: \[\int_{0}^\infty x\frac{1}{\sigma}\exp(-\frac{x}{\sigma})\] \[ = \sigma\]
Now plugging this back into \(**\) above we have the assymptotic variance is given by:
\[\frac{\sigma^2}{n}\]
35
Let \(U_1, U_2, \cdots U_{1029}\) be independent Uniformly distributed random variables. Let \(X_1\) be the number of \(U_i\) less than .331, \(X_2\) be the number of between .331 and .820 and \(X_3\) equal the number greater than .820. Why is the joint distribution of \(X_1, X_2, X_3\) multinomial with probabilities .331, .489, and .180 and n=1029?
Solution By definition of the uniform distribution there is equal probability for any value between zero and one. And so it follows that the probability of observing a value less than .331 is .331 and so on.
47.
The Pareto distribution has been used in economics as a model for a density function with a slowly decaying tail: \[f(x|x_0, \theta) = \theta x_0^\theta x^{-\theta-1}, \;\; x\geq x_0, \theta > 1\]
Assume that \(x_0 > 0\) is given and that \(X_1, X_2, ..., X_n\) is an i.i.d sample.
- Find the method of moments estimate of \(\theta\)
- Find the mle of \(\theta\)
- Find the asymptotic variance of the mle
Soluton
- First we find the theoretical moment of the pareto distribution call it \(\mu_1\)
\[E(X) = \int_{x_0}^\infty x\theta x_0^{\theta} x^{-\theta - 1}dx\] \[ = \theta x_0^\theta \Big[ \frac{x^{-\theta+1}}{-\theta + 1} \Big] \Big|_{x_0}^{\infty} = \frac{x_0^\theta}{\theta - 1}\]
\[\Rightarrow \theta = \frac{\mu}{\mu - x_0}\] and so the mm estimate becomes: \[\hat{\theta} = \bar{X}/(\bar{X} - x_0)\]
The log likelihood function is given by: \[l(\theta) = n\log(\theta) + n\theta \log(x_0) - (\theta + 1)\sum\log(x_i)\] and so taking the first derivative, setting to zero and solving we get: \[l'(\theta) = \frac{n}{\theta} + n\log(x_0) - \sum\log(x_i)\] \[l'(\theta) = \frac{n}{\theta} + n\log(x_0) - \sum\log(x_i) = 0\] \[\Rightarrow \hat{\theta} = \frac{n}{\sum\log(X_i) - n\log(x_0)}\]
The asymptotic variance is given by: \[1/[nI(\theta)]\]
So we find the Fisher Information: \[I(\theta) = -E[l''(\theta)] = E[\frac{\partial^2}{\partial\theta^2} (n\log(\theta) + n\theta \log(x_0) - (\theta + 1)\sum\log(x_i))]\] \[ = \frac{1}{\sigma^2}\] and so the asymptotic variance is: \[\sigma^2/n\]
46.
### Problem 46: Bowhead whale.
# a.
whale = scan("~/Downloads/whale.txt") # read data
whale.data = whale
br = seq(0,5,0.25) # create breaks for the histogram
X11()
hist(whale,breaks=br,main='whale data',probability = T,plot=T) # create histogram with breaks br
# since thew whale data is continuous, really should use a density estimator
#X11()
#plot(density(whale, n=50, window="g"),type="l",xlab="x",ylab="density")
########################################################################################
# matrices for output
out.mme = matrix(0,2,5)
out.mle = matrix(0,2,7)
########################################################################################
# b. MME
n = length(whale) # sample size
x.bar = mean(whale) # sample mean
sigmasq.hat = ((n-1)/n)*var(whale) # (1/n)sum(x_i - x.bar)^2
lambda.hat.mme = x.bar/sigmasq.hat
lambda.hat.mme## [1] 1.329121alpha.hat.mme = x.bar^2/sigmasq.hat
alpha.hat.mme## [1] 0.8002543out.mme[1,1] = alpha.hat.mme
out.mme[2,1] = lambda.hat.mme
# c. MLE
# solve the nonlinear equation
alpha.d = function(a){
n*log(a) - n*(log(mean(whale))) + sum(log(whale)) - n*digamma(a)
}
a.min = 0.1
a.max = 5
alpha.hat.mle = uniroot(alpha.d, c(a.min,a.max))
alpha.hat.mle = alpha.hat.mle$root
alpha.hat.mle## [1] 1.611355lambda.hat.mle = alpha.hat.mle / mean(whale)
lambda.hat.mle## [1] 2.676257out.mle[1,1] = alpha.hat.mle
out.mle[2,1] = lambda.hat.mle
# d.
# plot the fitted density using mme
par(new=T) # plot fitted curve on the histogram
x = seq(0.01,5,0.01)
y = dgamma(x,alpha.hat.mme,lambda.hat.mme)
y.max = 1.8 # approximate maximum of the histogram, used to match y-axis scale.
plot(x,y,type='l',ylim=c(0,y.max))
# plot the fitted density using mle
par(new=T) # plot fitted curve on the histogram
x = seq(0.01,5,0.01)
y = dgamma(x,alpha.hat.mle,lambda.hat.mle)
y.max = 1.8 # approximate maximum of the histogram, used to match y-axis scale.
plot(x,y,type='l',lty=4,ylim=c(0,y.max))# e. Bootstrap estimate of the s.e. and sampling distribution of the mme's
B = 1000
y = matrix(0,B,n)
lambda.star.mme = numeric(B)
alpha.star.mme = numeric(B)
for(i in 1:B){
y[i,] = rgamma(n, alpha.hat.mme,lambda.hat.mme)
x.bar = mean(y[i,])
sigmasq.hat = ((n-1)/n)*var(y[i,])
lambda.star.mme[i] = x.bar/sigmasq.hat
alpha.star.mme[i] = x.bar^2/sigmasq.hat
}
br = seq(0,5,0.10) # create breaks for the histogram
X11()
hist(alpha.star.mme,breaks=br,probability = T,plot=T) alpha.star.mean = mean(alpha.star.mme)
alpha.star.mean## [1] 0.8232669alpha.star.sd = sqrt(var(alpha.star.mme))
alpha.star.sd## [1] 0.1098088out.mme[1,2] = alpha.star.mean
out.mme[1,3] = alpha.star.sd
alpha.LCL = quantile(alpha.star.mme, 0.025)
alpha.LCL## 2.5%
## 0.6244449alpha.UCL = quantile(alpha.star.mme, 0.975)
alpha.UCL## 97.5%
## 1.035594out.mme[1,4] = alpha.LCL
out.mme[1,5] = alpha.UCL
X11()
hist(lambda.star.mme,breaks=br,probability = T,plot=T)lambda.star.mean = mean(lambda.star.mme)
lambda.star.mean## [1] 1.378505lambda.star.sd = sqrt(var(lambda.star.mme))
lambda.star.sd## [1] 0.2087042out.mme[2,2] = lambda.star.mean
out.mme[2,3] = lambda.star.sd
lambda.LCL = quantile(lambda.star.mme, 0.025)
lambda.LCL## 2.5%
## 0.9966149lambda.UCL = quantile(lambda.star.mme, 0.975)
lambda.UCL## 97.5%
## 1.810496out.mme[2,4] = lambda.LCL
out.mme[2,5] = lambda.UCL
# f. Bootstrap estimate of the s.e. and sampling distribution of the mle's
B = 1000
y = matrix(0,B,n)
lambda.star.mle = numeric(B)
alpha.star.mle = numeric(B)
a.min = 0.5
a.max = 5
for(i in 1:B){
y[i,] = rgamma(n, alpha.hat.mle,lambda.hat.mle)
whale = y[i,]
alpha.star.mle[i] = uniroot(alpha.d, c(a.min,a.max))$root # alpha
lambda.star.mle[i] = alpha.star.mle[i] / mean(whale) # lambda
}
br = seq(0,5,0.05) # create breaks for the histogram
X11()
hist(alpha.star.mle,breaks=br,probability = T,plot=T) alpha.star.mean = mean(alpha.star.mle)
alpha.star.mean## [1] 1.631555alpha.star.sd = sqrt(var(alpha.star.mle))
alpha.star.sd## [1] 0.148771out.mle[1,2] = alpha.star.mean
out.mle[1,3] = alpha.star.sd
alpha.LCL = quantile(alpha.star.mle, 0.025)
alpha.LCL## 2.5%
## 1.367639alpha.UCL = quantile(alpha.star.mle, 0.975)
alpha.UCL## 97.5%
## 1.94607out.mle[1,4] = alpha.LCL
out.mle[1,5] = alpha.UCL
X11()
hist(lambda.star.mle,breaks=br,probability = T,plot=T)lambda.star.mean = mean(lambda.star.mle)
lambda.star.mean## [1] 2.721388lambda.star.sd = sqrt(var(lambda.star.mle))
lambda.star.sd## [1] 0.2894131out.mle[2,2] = lambda.star.mean
out.mle[2,3] = lambda.star.sd
lambda.LCL = quantile(lambda.star.mle, 0.025)
lambda.LCL## 2.5%
## 2.217553lambda.UCL = quantile(lambda.star.mle, 0.975)
lambda.UCL## 97.5%
## 3.39229out.mle[2,4] = lambda.LCL
out.mle[2,5] = lambda.UCL
# g. approximate 95% CI for the mle's
# alpha
delta.low = quantile(alpha.star.mle,0.025) - alpha.star.mean
delta.up = quantile(alpha.star.mle,0.975) - alpha.star.mean
mle.LCL = alpha.hat.mle - delta.up
mle.LCL## 97.5%
## 1.296841mle.UCL = alpha.hat.mle - delta.low
mle.UCL## 2.5%
## 1.875272out.mle[1,6] = mle.LCL
out.mle[1,7] = mle.UCL
# lambda
delta.low = quantile(lambda.star.mle,0.025) - lambda.star.mean
delta.up = quantile(lambda.star.mle,0.975) - lambda.star.mean
mle.LCL = lambda.hat.mle - delta.up
mle.LCL## 97.5%
## 2.005355mle.UCL = lambda.hat.mle - delta.low
mle.UCL## 2.5%
## 3.180093out.mle[2,6] = mle.LCL
out.mle[2,7] = mle.UCL
# Summary
dimnames(out.mme) = list(c("alpha","lambda"),c("MME", "Bootstrap.mean", "Bootstrap.sd", "quant.LCL", "quant.UCL"))
dimnames(out.mle) = list(c("alpha","lambda"),c("MLE", "Bootstrap.mean", "Bootstrap.sd", "quant.LCL", "quant.UCL", "approx.LCL", "approx.UCL"))
out.mme## MME Bootstrap.mean Bootstrap.sd quant.LCL quant.UCL
## alpha 0.8002543 0.8232669 0.1098088 0.6244449 1.035594
## lambda 1.3291210 1.3785047 0.2087042 0.9966149 1.810496out.mle## MLE Bootstrap.mean Bootstrap.sd quant.LCL quant.UCL approx.LCL
## alpha 1.611355 1.631555 0.1487710 1.367639 1.94607 1.296841
## lambda 2.676257 2.721388 0.2894131 2.217553 3.39229 2.005355
## approx.UCL
## alpha 1.875272
## lambda 3.18009350.
Solution
Find the method of moments estimator: Finding the theoretical moment we have: \[E(X) = \int_0^\infty \frac{x^2}{\theta^2}e^{\frac{-x^2}{2\theta^2}}\] and so we have that the estimate is given by: \[\hat{\theta} = \sqrt{\frac{2}{\pi}\bar{X}}\]
To find the mle we write out the log ligekilihood function, which is given by: \[l(\theta) = \sum[\ln(x_i) - 2\ln(\theta) - \frac{x_i^2}{2\theta^2}]\] \[ = \sqrt{\frac{\pi\theta^2}{2}}\] now taking the first derivative, setting to 0 and solving we have: \[l'(\theta) = -2n/\theta + \sum x_i^2/\theta^3\] \[\Rightarrow \hat{\theta} = \sqrt{\frac{1}{2n}\sum X_i^2}\]
We simply use what was formulated above and find that the asymptotic variance is: \[\frac{\theta^2}{4n}\]
52.
- Find the asymptotic variance of the mle
Solution
- We know from previous results that: \[l''(\theta) = -n/(1+\theta)^2\] and so it follows that the asymptotic variance is given by: \[(1 + \theta)^2/n\]
55.
Solution
(a/b) To find the mle we write out the log-likehood function given by:
\[\ln(n!) - \sum\ln(X_i!) + X_1\ln .25(2 + \theta) + X_2\ln .25(1-\theta) + X_3\ln .25(1- \theta) + X_4\ln .25\theta\] Now taking the derivative of this equating to zero and solving yields: \[\theta^2(-X_1 - X_2 - X_3 - X_4) + \theta(X_1 - 2X_2 - 2X_3 - X_4) + 2X_4 = 0\] Using the provided data we get: \[0.035712\] Moreover using the formula we find the asymptotic variance to be: \[1/110076.9 \]
Using part (a) we can form the 95% cofidence interval to get: \[(0.02980475,0.04161986)\]
p1 <- 0.25*(2+theta_hat)
p2 <- 0.25*(1-theta_hat)
p3 <- 0.25*(1-theta_hat)
p4 <- 0.25*theta_hat
numberOfReps <- 1000
sim <- rmultinom(numberOfReps, size = n, prob = c(p1,p2,p3,p4))
est <- c()
for(i in 1:numberOfReps){
X1 <- sim[1,i]
X2 <- sim[2,i]
X3 <- sim[3,i]
X4 <- sim[4,i]
cT2 <- -X1-X2-X3-X4
cT <-X1-2*X2-2*X3-X4
C <- 2*X4
roots <- as.numeric(polyroot(c(C, cT, cT2)))
est[i] <- roots[roots >= 0]
}
sd(est)which results in $ 0.005824132$
58.
On attached handwritten paper.