The University of Sydney
School of Mathematics and Statistics
Computer Exercise Week 10
STAT3023: Statistical Inference Semester 2, 2021
Lecturers: Rachel Wang and Michael Stewart
Interval estimation of a Binomial success probability
Suppose we have a single B(m, θ) observation X for m known but the success probability parameter
θ ∈ Θ = (0, 1) unknown. Consider the decision problem where the decision space is D = Θ and the loss
function L(d|θ) = 1 {|d− θ| > C}. Then a “decision” d(X) may be viewed as an estimator of θ and the
risk is precisely the probability that the resultant interval estimator d(X)±C does not cover the true θ,
i.e. the noncoverage probability of the interval. In other words, we are interested in forming a fixed-width
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).
We shall consider two different procedures:
• a simple procedure centred at the observed proportion θˆML = Xm ;
• a Bayes procedure using the “flat” prior w(θ) ≡ 1; note that in this case this is a “proper prior”,
in fact the beta(1, 1) density, an example of a “conjugate prior” for this problem.
For simplicity, we shall take m = 10 and suppose we wish to form an interval of width 0.2, so C = 0.1.
Since there are only 11 different possible observations, we can exhaustively list all possible intervals
using both procedures. It is then possible to compute “exactly” the risk functions of the two procedures.
Note also that when X takes the (extreme possible) values 0 and 10, we have a modified rule for
constructing the interval, to ensure it remains within (0,1).
1. The simple interval
For observed values x = 1, 2, . . . , 9, the simple interval is
[
x−1
10 ,
x+1
10
]
. But for x = 0 the interval is
the same as for x = 1; for x = 10 the interval is the same as for x = 9.
(a) Define x=0:10 and also define two vectors simple.lower and simple.upper giving the lower
and upper endpoints of the simple interval.
(b) Define now th=(1:1000)/1001. The following code computes the risk of the simple interval
for each value in th:
th=(1:1000)/1001
L=length(th)
risk=0
for (i in 1:L){
probs=dbinom(x,10,th[i])
lower.miss=th[i]upper.miss=th[i]>simple.upper
risk[i]=sum(probs[lower.miss])+sum(probs[upper.miss])
}
Plot risk against th (as a red line plot).
2. The Bayes interval
The product of the likelihood and the weight function/prior is(
m
X
)
θX(1− θ)m−X = 1
m+ 1︸ ︷︷ ︸
const.
θ(X+1)−1(1− θ)(m−X+1)−1
B(X + 1,m−X + 1)︸ ︷︷ ︸
p(θ|X)
where the B(X + 1,m −X + 1) = X!(m−X)!(m+1)! in the denominator here denotes the beta function;
the posterior is thus the beta(X + 1,m−X + 1) density.
Copyright© 2021 The University of Sydney 1
The interval is determined by finding the level set of the posterior density of width 0.2. For X = 0,
the posterior is just the decreasing function
(m+ 1)(1− θ)m
and so the level set is simply [0, 0.2]. Similarly for X = 10, the interval is [0.8, 1].
For all other values X = 1, 2, . . . , 9, the posterior density increases and then decreases and we need
to use the computer (specifically, the R function uniroot()) to find the level set; that is, we need
to find the value d such that
p(d− C|X) = p(d+ C|X) ;
then the desired interval is d± C.
(a) First define an objective function for which to find the root
beta.obj=function(d,X,m=10,C=0.1){
dbeta(d-C,X+1,m-X+1)-dbeta(d+C,X+1,m-X+1)
}
(b) Try solving this manually for X = 1:
uniroot(beta.obj,lower=0.1,upper=0.9,X=1)
(c) Verify that we have indeed found the right value:
• extract the root;
• form the interval by defining lower.endpt and upper.endpt;
• plot the posterior (beta(2,10)) density (try using curve());
• add vertical lines (using lty=2) at lower.endpt and upper.endpt;
• add a horizontal line (using lty=3) showing the level set (i.e. at the height p(θ|X)
where θ = d± C and X = 1).
(d) Now iterate through the values X = 1, . . . , 9 to find the lower and upper endpoints of
the Bayes intervals, saving these in vectors bayes.lower and bayes.upper. Complete the
vectors with appropriate endpoints for X = 0 and X = 10.
(e) Adapt the code from question 1 part (b) above to compute the risk of the Bayes interval.
Recreate your plot from 1(b) and add to it the Bayes interval risk in blue. Add an informative
heading and legend.
3. For which values of th
• does the simple interval do better;
• does the bayes interval do better;
• do the two intervals have similar performance?
2

欢迎咨询51作业君