This is an exam. You can use any written material, either in books or online but you can not discuss it with ANYONE.
Do as much as possible analytically, otherwise use numerical method or simulation.

Below is a sample from a population with density \(f(x\vert a)=a/x^{a+1};x>1, a>1\)

1 1.02 1.03 1.03 1.03 1.03 1.03 1.04 1.05 1.05 1.05 1.06 1.06 1.07 1.07 1.08 1.09 1.09 1.09 1.1 1.11 1.11 1.12 1.13 1.13 1.14 1.15 1.15 1.16 1.16 1.16 1.17 1.2 1.2 1.21 1.21 1.22 1.23 1.23 1.23 1.24 1.24 1.24 1.24 1.24 1.25 1.25 1.26 1.28 1.34 1.35 1.41 1.43 1.43 1.43 1.44 1.44 1.44 1.44 1.44 1.45 1.45 1.46 1.47 1.48 1.48 1.51 1.51 1.53 1.55 1.56 1.58 1.58 1.61 1.65 1.66 1.69 1.71 1.73 1.75 1.77 1.78 1.8 1.84 1.87 1.88 1.93 1.97 2.02 2.12 2.26 2.3 2.52 2.62 2.7 2.73 3.1 3.74 3.82 4.29

You can generate your own samples with

rmid=function(n,a) runif(n)^(-1/a)
  1. Find the method of moments and the maximum likelihood estimators of a.

\[ \begin{aligned} &E[X] =\int_1^\infty xa/x^{a+1}dx=\int_1^\infty ax^{-a}dx= \\ &\frac{a}{-a+1}x^{-a+1}\vert_1^\infty = \frac{a}{a-1}\\ &\frac{a}{a-1} =\bar{x} \\ &\hat{a}=\frac{\bar{x}}{\bar{x}-1} \end{aligned} \]

\[ \begin{aligned} &f(\pmb{x}\vert a) =a^n/\left(\prod_{i=1}^n x_i\right)^{a+1} \\ &l(a\vert\pmb{x}) =n\log a-(a+1)\sum_{i=1}^n\log x_i \\ &dl(a\vert\pmb{x})/da =n/a-\sum_{i=1}^n\log x_i=0 \\ &\hat{a}=1/\overline{\log x} \end{aligned} \]

round(c(mean(dta)/(mean(dta)-1), 1/mean(log(dta))), 3)
## [1] 2.931 2.764
  1. Find the bias of the mle.

\[ \begin{aligned} &F(x\vert a)=\int_1^x a/t^{a+1}dt=-t^{-a}\vert_1^x = 1-1/x^a;a>1,x>1\\ &P(\log X<x) = P(X<e^x)=1-1/(e^x)^a=1-e^{-ax};x>1\\ &\log X\sim Exp(a) \\ &T=\sum_{i=1}^n \log X_i \sim Gamma(n, 1/a) \\ &E[\hat{a}]=E[n/T]=\int_{1}^{\infty}\frac{n}{x} \frac{a^n}{\Gamma(n)}x^{n-1}e^{-ax}dx=\\ &an\Gamma(n-1)/\Gamma(n)\int_{1}^{\infty} \frac{a^{n-1}}{\Gamma(n-1)}x^{(n-1)-1}e^{-ax}dx=an/(n-1) \end{aligned} \] Now

\[bias(\hat{a})=E[\hat{a}]-a=an/(n-1)-a=a/(n-1)\]

Here is a little simulation that shows the same:

B=1e4
a=2;n=3
x=matrix(rmid(n*B, 2), B, n)
ahat=1/apply(log(x), 1, mean)
round(c(a, mean(ahat), a*n/(n-1), a/(n-1)), 3)
## [1] 2.000 3.052 3.000 1.000
  1. Is the mle a consistent estimator of a?

\[ \begin{aligned} &P\left(\vert\hat{a}-a\vert<\epsilon\right) = \\ &P\left(\vert n/T-a\vert<\epsilon\right) = \\ &P\left(a-\epsilon< n/T<a+\epsilon\right) = \\ &P\left(n/(a+\epsilon)< T<n/(a-\epsilon)\right) = \\ &\int_{n/(a+\epsilon)}^{n/(a-\epsilon)} \frac{a^n}{\Gamma(n)}x^{n-1}e^{-ax}dx =\\ &\int_{n/(a+\epsilon)}^{n/(a-\epsilon)} \frac{1}{\Gamma(n)}(ax)^{n-1}e^{-ax}(adx) = \\ &\int_{an/(a+\epsilon)}^{an/(a-\epsilon)} \frac{1}{\Gamma(n)}t^{n-1}e^{-t}dt \end{aligned} \] there does not seem to be an easy way to find this integral, so we might use numerical method:

n=10000;a=2;eps=0.1
I=a*n/c(a+eps, a-eps)
diff(pgamma(I, n, 1))
## [1] 0.9999992

and so it appears the mle is consistent.

  1. Using the p value approach do the likelihood ratio test at the 5% level to test

\[H_0:a=2\text{ vs. }H_a:a\ne 2\] based on the large sample approximation of the likelihood ratio test

\[ \begin{aligned} &2(l(\hat{a})-l(a_0)) = \\ &2\left(n\log \hat{a}-(\hat{a}+1)\sum_{i=1}^n\log x_i-n\log a_0+(a_0+1)\sum_{i=1}^n\log x_i\right) = \\ &2\left(n\log (n/T)-(n/T+1)T-n\log a_0+(a_0+1)T\right) = \\ &2n\log [n/(a_0T)]- 2(n-Ta_0) \end{aligned} \]

n=length(dta);a0=2
TS=sum(log(dta))
ahat=n/TS
lrt=2*n*log(n/TS/a0)-2*(n-a0*TS)
round(c(lrt, 1-pchisq(lrt, 1)), 4)
## [1] 9.4279 0.0021

and so we fail to reject the null hypothesis.

  1. Do the likelihood ratio test at the 5% level to test

\[H_0:a=2\text{ vs. }H_a:a\ne 2\] without using the large sample approximation of the likelihood ratio test.

let’s draw the curve of the likelihood ratio test statistic. Note that if a=3 and n=100 \(E[T]=n/a=100/3\) and \(sd[T]=\sqrt{n}/a=10/3\).

lrt=function(x) 2*n*log(n/x/a0)-2*(n-a0*x)
curve(lrt, 20, 50)

so clearly “lrt is small” is equivalent to “T” is neither small nor large. So

\[ \begin{aligned} &\alpha =P\left(|T|>cr\right) \\ &\alpha/2 =P\left(T<cr_1\right) \\ &\alpha/2 =P\left(T>cr_2\right) \\ &cr_1=qgamma(\alpha/2,n, a_0)\\ &cr_2=qgamma(1-\alpha/2,n, a_0) \end{aligned} \]

n=length(dta);alpha=0.05
qgamma(c(alpha/2,1-alpha/2), n, a0)
## [1] 40.68200 60.26447
TS
## [1] 36.17702

and so again we reject the null hypothesis.

  1. Find the Bayesian point estimator of a using the posterior mean and the prior \(\pi(a)=1/a;a>1\)

\[ \begin{aligned} &f(\pmb{x}, a) = a^n/\left(\prod_{i=1}^n x_i\right)^{a+1}\times 1/a = a^{n-1} e^{-(a+1)T}\\ &m(T) =\int_{1}^\infty a^{n-1} e^{-(a+1)T}da = \\ &e^{-T}\frac{\Gamma(n)}{T^n}\int_{1}^\infty \frac{1}{\Gamma(n)}(Ta)^{n-1} e^{-aT}(Tda) = \\ &e^{-T}\frac{\Gamma(n)}{T^n}\int_{T}^\infty \frac{1}{\Gamma(n)}t^{n-1} e^{-t}dt=\\ &e^{-T}\frac{\Gamma(n)}{T^n}\left(1-pgamma(T, n, 1)\right)\\ &f(a\vert T) = \frac{a^{n-1} e^{-(a+1)T}}{e^{-T}\frac{\Gamma(n)}{T^n}\left(1-pgamma(T, n, 1)\right)} = \\ &\frac{T^n}{\Gamma{(n)}}a^{n-1}a^{-aT}/\left(1-pgamma(T, n, 1)\right) \end{aligned} \] if \(G\sim\Gamma(n,T)\), then the posterior density is the density of \(G|G>1\).

Finally

\[ \begin{aligned} &E[a\vert T] = \int_1^\infty a \frac{T^n}{\Gamma{(n)}}a^{n-1}a^{-aT}/\left(1-pgamma(T, n, 1)\right) da =\\ &\frac{\Gamma(n+1)}{T\Gamma(n)}\int_1^\infty \frac{T^{n+1}}{\Gamma{(n+1)}}a^{(n+1)-1}a^{-aT}/\left(1-pgamma(T, n, 1)\right) da = = \\ &\frac{n}{T}\frac{1-pgamma(T, n+1, 1)}{1-pgamma(T, n, 1)} \end{aligned} \]

n=length(dta)
TS=sum(log(dta))
(n/TS)*(1-pgamma(TS, n+1, 1))/(1-pgamma(TS, n, 1))
## [1] 2.764185