Solution to Homework 4

Say \(X_1,..,X_n \sim Exp(\lambda)\). Let \(T = 1/\bar{x}\)

1. find \(Z\), the Bayesian estimator of \(\lambda\) when \(\lambda\sim Exp(1)\) and using the mean of the posterior distribution.

2. which of these estimators is unbiased?

\(X_1,..,X_n\) iid Exp( \(\lambda\) ), so \(\sum X_i\sim Gamma(n,1/ \lambda)\). Therefore

and so \(T\) is (slightly) biased.

Now

and this integral can not be found analytically. Here is a numerical example that shows this estimator is not unbiased:

f <- function(x, l=1, n=3) {
  (n+1)*l^n/factorial(n-1)/(x+1)*x^(n-1)*exp(-l*x)
}
integrate(f, 0, Inf)$value

## [1] 1.192695

3. which of these estimators is a sufficient statistic?

both of them are sufficient statistics.

4. find a function \(g(T)\) which is an ancillary statistic of \(\lambda\) .

`

now the right hand side does not depend on \(\lambda\) , and so \(g(T) = T/ \lambda\) is an ancillary statistic for \(\lambda\) .

Problem 5 : Is \(T\) a consistent estimator of \(\lambda\) ?

6. consider the estimators \(T\) and \(V=\log (2)/\text{median}(x)\). Write a routine in R that uses simulation to draw a curve for eff(V|T) either as a function of n (with \(\lambda\) fixed) or as a function of \(\lambda\) (with n fixed). (You can assume that T and V are unbiased)

hw4p6 <- function (problem, n, lambda, B=1e4, npoints=50) {
  if(length(n)==1) {
  x <- seq(lambda[1], lambda[2], length=npoints)
  y <- rep(0,npoints)
  for(i in 1:npoints) {
    A <- matrix(rexp(n*B, x[i]), ncol=n)
    U <- (n-1)/apply(A, 1, sum)
    V <- log(2)/apply(A, 1, median)
    y[i] <- var(V)/var(U)
  } 
 }
 else {
  n <- n[1]:n[2]
  y <- 0*n
  for(i in 1:length(n)) {
    x <- matrix(rexp(n[i]*B,lambda),ncol=n[i])
    U <- (n[i]-1)/apply(x, 1, sum)
    V <- log(2)/apply(x, 1, median)
    y[i] <- var(V)/var(U)
  }
  x <- n
 }
 plt <- ggplot(data.frame(x=x, y=y), aes(x, y)) +
   ylab("Relative Efficiency")
  if(length(n)==1) 
    plt <- plt + geom_line() + xlab(expression(lambda))
 else
   plt <- plt + geom_point() + xlab("n")
 print(plt)
 
}

hw4p6(n=c(10, 50), lambda=1.5)

hw4p6(n=20, lambda=c(1.5, 10))