STAT3023/3923/4023: Statistical Inference Semester 2, 2022

Suppose X =  (X1 , . . . , Xn ) consists of iid Poisson(θ) random variables with unknown mean θ  e R. Consider the decision problem of estimating θ with squared-error loss,

L(dlθ) = (d _ θ)2

so that risk of an estimator d(X) is the mean-squared error :

R(θld) = Eθ { [d(X) _ θ]2 } .

Consider two estimators:

_  The maximum likelihood estimator dmle (x) =  (the sample mean).

_  The Bayes procedure dconj (x) using a unit-mean exponential prior (note that this is an example

of a conjugate prior). The product of the prior and the likelihood is

n    e −θ θ Xi

i=1      Xi !

where T =      Xi  = n is the sample total.  Since (as a function of θ) this is proportional to a gamma density with shape T + 1 and rate n + 1, the posterior mean is

shape      T + 1

rate       n + 1 .

We are going to approximate rescaled risk via simulation and plot it, that is we are going to _  define a vector of θ values;

_  for each θ value, we will

– generate B samples of size n from that Poisson(θ) distribution, and obtain the value of the estimators dmle (x) =  and dconj (x) for each;

– obtain the average squared error for each estimator;

_ we will then plot n times the average squared error against θ for both estimators. We are going to do this for two different sample sizes: n=5 and n=50.

1.  Define a vector of θ values using th=1:100/10.  Define also L .th=length(th), n=5 and define empty vectors nMSE .mle=0 and nMSE .conj=0. Finally define B=10000 as the number of simulation iterations for each θ value (if you computer is very slow you can try making B smaller, say B=1000).

2.  Write a “double loop” . At the j-th iteration of the outer loop,

_  initialise mle=0, conj=0;

_  simulate B samples of size n from the Poisson distribution with mean th[j], obtaining values

of each estimator, saving them in the vectors mle and conj respectively;

_  compute n times the average mean-squared error between the estimates in mle and th[j]

and save it in nMSE .mle[j]; do the same for conj and nMSE .conj[j].