The University of Sydney
School of Mathematics and Statistics
Computer Exercise Week 11
STAT3023: Statistical Inference Semester 2, 2021
Lecturers: Neville Weber and Michael Stewart
Prepare your computer report and submit to the appropriate Canvas portal by 11:59pm
Sunday 31st October. Please include in your report all the code, plots and any comments
required by the questions. The permitted format of the file that you upload will be restricted
to PDF, HTML or Word, and should be created in a reproducible fashion (e.g. compiled from
an Rmarkdown file). Please do not have your name visible in your report.
Interval estimation of an exponential rate parameter
Suppose we have a small sample of iid observations from an exponential distribution with unknown
rate θ ∈ Θ = (0,∞). Consider the formal decision problem where the decision space is D = Θ = (0,∞)
and the loss function is L(d|θ) = 1 {|d− θ| > C} for some fixed, known C > 0. In other words, we are
interested in forming a fixed-width (of 2C) interval estimate of θ (rather than a “confidence interval”;
the main difference is that we are fixing the width, not the coverage probability of the procedure). The
risk is the non-coverage probability of this interval; the estimator is the midpoint of the interval.
We are interested in the risk functions of two “decisions”.
1. An interval based on the maximum likelihood estimator θ̂ML = 1/X̄: a “näıve” interval would
simply be θ̂ML ± C, although this would not make complete sense if θ̂ML < C since then part of the interval would be below zero; we thus modify the interval by defining the centre as d1(X) = max ( θ̂ML, C ) = { θ̂ML if θ̂ML ≥ C C otherwise. Then the interval estimate is d1(X)± C = {( θ̂ML − C, θ̂ML + C ) if θ̂ML ≥ C (0, 2C) otherwise. 2. A Bayes procedure for this decision problem using the weight function w(θ) ≡ 1 (the “flat prior”). Since the likelihood is fθ(X) = n∏ i=1 [ θe−θXi ] = θne−Tθ , where T = n∑ i=1 Xi, the product of the likelihood and prior/weight function is a multiple of a gamma density with shape n+ 1 and rate T ; that is, the posterior density is p(θ|X) = θ(n+1)−1e−TθTn+1 Γ(n+ 1) . The Bayes estimator d2(X) is the midpoint of the level set of this density of width 2C; that is, d2 = d2(X) satisfies p(d2 − C|X) = p(d2 + C|X) . (1) Equivalently, (d2 − C) n e−T (d2−C) = (d2 + C) n e−T (d2+C) e2TC = ( d2 + C d2 − C )n e 2TC n = d2 + C d2 − C (d2 − C)e2X̄C = d2 + C d2 = C ( e2X̄C + 1 e2X̄C − 1 ) . Copyright© 2021 The University of Sydney 1 The endpoints of the interval are given by d2 − C = 2C e2X̄C − 1 , and d2 + C = 2Ce2X̄C e2X̄C − 1 . For simplicity we shall take n = 4 and C = 1.5. We shall make use of the fact that X̄ has a gamma distribution with shape n and rate nθ. We shall firstly approximate the risk functions via simulation, then compute the exact risk functions. 1. (a) Write a function mle.int() that computes the mle-based interval. It should look like mle.int=function(x,C) { ... d1=... c(d1-C,d1+C) } where, of course, the ... parts are replaced with proper R code! (b) Test it out by generating a test sample: n=4 C=1.5 th0=2 x=rexp(n,rate=th0) mle.int(x,C) 2. (a) Now write a function bayes.int() which computes the Bayes interval based on the flat prior. It should look like bayes.int=function(x,C) { ... d2=... c(d2-C,d2+C) } (b) Try it out on your test sample to compare the two intervals. (c) Visualise the interval by • plotting the posterior density (using curve()); • adding vertical lines (using lty=2) indicating the interval endpoints; • adding a horizontal line (using lty=3) at the (common) height of the posterior density at the two interval endpoints; (see week 10 2(c) for a similar plot). Hint: define m=mean(x), use n and m to define the shape and rate of the posterior (gamma) density, and use curve(...,from=0,to=3/m). 3. Define th=(1:500)/50, L=length(th) and B=1000. Also define noncoverage.mle=noncoverage.bayes=0 Perform a double loop: at the i-th iteration of the outer loop • define matrices mle.mat=matrix(0,B,2) and bayes.mat=matrix(0,B,2); • perform the inner loop: at the j-th iteration of the inner loop – generate a pseudo-random sample x of size n=4 from an exponential distribution with rate th[i]; 2 – obtain mle-based and Bayes intervals (with C=1.5), saving them in the j-th row of mle.mat and bayes.mat respectively; • save in the i-th element of noncoverage.mle the number of times the mle-based interval did not cover th[i]; similarly for the Bayes interval in bayes.mat; Convert the counts in the noncoverage vectors to proportions. Plot (as lines) these proportions against th (red for the mle-based, blue for Bayes). Add an informative heading and legend. 4. The precise form of the risk function for the mle-based interval is different for 0 < θ < 2C and θ ≥ 2C, because for 0 < θ < 2C, the interval can only miss if the mle θ̂ML is too large, whereas for θ ≥ 2C, the interval can miss on both the high side and the low side. Therefore, for 0 < θ < 2C, R(θ|d1) = Pθ { θ < θ̂ML − C } = Pθ { X̄ < 1 θ + C } and for θ ≥ 2C, R(θ|d1) = Pθ { θ < θ̂ML − C } + Pθ { θ > θ̂ML + C
}
= Pθ
{
X̄ <
1
θ + C
}
+ Pθ
{
X̄ >
1
θ − C
}
.
The first probability here is of overestimating, and applies to all values of θ. The second is the
probability of underestimating and only applies to θ ≥ 2C.
Since X̄ has a gamma distribution with shape n and rate nθ, we can compute this risk exactly.
Compute this risk function for each value in th, save it in a vector called m.risk and plot it against
th. Your code should look something like this:
risk.overest = … # this applies to all values in th
big = (th>=(2*C)) # logical indicator of th values >= 2*C
risk.underest = 0*risk.overest # start with a vector of zeroes
risk.underest[big] = … # something involving th[big]
m.risk = risk.overest + risk.underest
5. The risk of the Bayes interval is given by
R(θ|d2) = Pθ(θ < d2 − C) + Pθ(θ > d2 + C) .
The first probability is
Pθ
(
θ <
2C
e2X̄C − 1
)
= Pθ
(
X̄ <
1
2C
log
(
1 +
2C
θ
))
.
Again, because d2 ≥ C (the Bayes interval’s lower endpoint never drops below zero!), the event in
the second probability can only happen for θ ≥ 2C. For such θ it is
Pθ
(
θ >
2Ce2X̄C
e2X̄C − 1
)
= Pθ
(
X̄ >
1
2C
log
(
θ
θ − 2C
))
.
Compute this risk function for each value in th (save it in a vector b.risk) and plot it (in
blue) against th (you will need to make use of the big indicator and compute risk.overest and
risk.underest separately again as you did in the previous question). Add a line plot of the risk
of the mle-based interval (in red). Add an informative heading and legend.
6. For which values of th
• does the mle-based interval do better;
• does the Bayes interval do better;
• do the two intervals have similar performance?
Can you describe and explain the behaviour of both curves near θ = 3?
3