Homework #4 for Stat 610

Due Thursday March 22, in class.

1. Recall the example from class.  X_1,...,X_n measure tree height theta,
   and we suppose X_i ~ Normal[theta, theta^2].  We found the MLE

   theta.hat <-  ( -bar X_n + sqrt{ 5 bar X_n^2 + 4 S_n^2  } )/2

   The following 10 heights were measured (ft) on a given tree:

   257.2, 48.8, 214.3, 17.0, 118.6, 178.9, 286.9, 186.1, 229.6, 113.2

   Compute the MLE and an approximate 95% confidence interval for theta.
   

2. Random variables X_1, X_2, ..., X_n are iid from a Negative Binomial
   distribution, with mean lambda>0 and shape alpha > 0.  Recall
   X_i may be constructed as a Poisson distributed variable with
   mean lambda*z_i, and with z_i the realization of 
   Z_i ~ Gamma( alpha, alpha ).  
   a. Derive expressions for E(X_i) and var(X_i).
   b. With bar X_n and S^2_n the sample mean and variance respectively,
      obtain a formula for a method of moments estimator of lambda
      and alpha.
   c. Write the likelihood, log likelihood, and score function for
      the parameters alpha and lambda. Show that the MLE of lambda
      is bar X_n.  Although there is no closed formula solution for 
      the MLE of alpha, a solution may be derived by Newton's method
      or some related method given a data set.

   d. Adapt the R source below to simulate B negative binomial 
      data sets, in order to compare the method of moments estimator
      of alpha with the MLE.

      #####################################################
      B <- 1000
      n <- 20
      alpha <- 3
      lambda <- 10
      alpha.mom <- numeric(B)
      alpha.mle <- numeric(B)

      for( b in 1:B )
       {
        z <- rgamma(n,shape=alpha,rate=alpha)
	x <- rpois(n,lambda=z*lambda)
	lambda.hat <- mean(x)
	alpha.mom[b] <- ** from part (b) above **
	alpha.mle[b] <- ** result of some optimization **
	print(b)
       }
      summary(alpha.mom)
      summary(alpha.mle)
      #####################################################


3. Consider a one-parameter exponential family in the natural
   parameterization. That is, X has density/mass function

    f(x;eta) = h(x) exp[ eta T(x) - A(eta) ]

   a. Establish an expression for var[ T(X) ] in terms of A(eta).
   b. A function g(eta) for which g''(eta) > 0 at all eta
      is necessarily convex.  Show that the negative loglikelihood for
      eta is a convex function.

4. a.  Show  that g(x) = x^p is convex when p >=1, for x >= 0.
   b.  Let x_1, x_2, ..., x_n denote positive real numbers.
       Show that 

         (1/n) sum_{i=1}^n x_i  >= ( prod_{i=1}^n x_i )^(1/n)

5. Consider data X and missing data Z, for which we have a joint
   probability model f(x,z;theta).  The log likelihood of interest is
   l(theta) = log f(x;theta), based on the marginal distribution of 
   X taken at the observed value x.  For an initial guess of the
   parameter, theta1, construct

   J(theta,theta1) = E_theta1[  log f( X, Z ; theta ) | X=x ]

   which is constructed during the E-step of the EM algorithm.
   Define

   h(theta,theta1) = J(theta,theta1) - J(theta1,theta1) + l(theta1)

   a. Show that

   h(theta,theta1) = E_theta1 { log [  f(X,Z;theta)/f(X,Z;theta1) ] | X=x } +
                     l(theta1)

   b. Use Jensen's inequality and facts about conditional distributions
   to show that

   h(theta,theta1) <= l(theta) for all theta.

   Since the M-step maximizes J, it also maximizes h.  Thus the
   M-step is an instance of optimization transfer, producing a 
   new point theta2 for which l(theta2) >= l(theta1).