# Statistics 2 Worksheet 9

Document Sample

```					                      Statistics 2: Worksheet 9
Jonathan Rougier, January 2010

This homework does not have to be handed in. Work through it carefully to improve
your understanding of the Bayesian approach to inference. There is a combination of
examples and questions; answers will be posted on Wed 20th January.

Remember the Dutch farmer from HW07?

A Dutch farmer buys a large box of tulip bulbs from a wholesaler. How-
ever, when he gets home, he cannot remember whether he bought the box
with (i) 60% red and 40% yellow bulbs, or (ii) 25% red and 75% yellow.
So he plants 12 bulbs. He decides he will select (i) if 8 or more come up
red, and (ii) otherwise.

This problem has the general form X ∼ Binomial(12, θ), where X is the number of
bulbs that come up red, and θ is the probability that an individual bulb will come up
red. Let’s suppose that xobs = 8: this is one way in which the Bayesian approach is
very diﬀerent from the Frequentist one: the only value of X that matters is xobs .

Discrete parameter

In general, θ ∈ [0, 1], but in our particular case, to start with, we will take θ ∈ Ω =
{0.25, 0.6}. What is the posterior probability distribution for θ? We use Bayes’s
Theorem, which tells us
q(θ) = m−1 L(θ) p(θ)

where p(θ) is the prior distribution for θ, and q(θ) is the posterior distribution. m−1
is the constant of proportionality, which can be found (for discrete Ω) as

m=           L(θ) p(θ)
θ∈Ω

We cannot ﬁnd the posterior distribution without specifying a prior distribution.
So let’s suppose, since we are given no other information, that the farmer thinks that
each outcome is equally likely, i.e.

0.5 θ ∈ {0.25, 0.6}
p(θ) =
0   otherwise.

For the likelihood we have, as usual,

obs               obs )
L(θ) ∝ f (xobs ; m, θ) ∝ θx (1 − θ)(m−x             = θ8 (1 − θ)4 .

1
Hence the posterior distribution has the form

m−1 × θ8 (1 − θ)4 × 0.5 θ ∈ {0.25, 0.6}
q(θ) =
0                       otherwise,

where
m=                  θ8 (1 − θ)4 × 0.5.
θ∈{0.25,0.6}

Here’s how we do this calculation in R. First, we deﬁne a vector of possible θ values,
which we will call Omega. Then we specify the prior probabilities for the components
of Omega.

> Omega <- c(0.25, 0.6)
> pOmega <- c(0.5, 0.5)

Now we specify the likelihood for each component of Omega:

> LOmega <- Omega^8 * (1 - Omega)^4

Finally, we compute the constant of proportionality, and from this derive the posterior
probabilities for Omega:

> m <- sum(LOmega * pOmega)
> qOmega <- LOmega * pOmega/m
> qOmega

[1] 0.01110365 0.98889635

So eight bulbs came up, and the result is that the posterior probability on θ = 0.25 is
1% and the posterior probability on θ = 0.6 is 99%. Our farmer is now almost certain
that he has bought the box with 60% red bulbs.

Q1. Repeat the analysis with a larger collection for Omega, say

> Omega <- seq(from = 0, to = 1, by = 0.1)

where, as before, the prior probability is the same for each possible value. If
you do this correctly, you should be able to draw the following ﬁgure:

2
Prior and posterior probabilities for theta

0.30
q
q     Prior
0.25
0.20   q     Posterior                        q
Probability

q
0.15

q
0.10

q       q     q    q    q         q    q    q    q     q   q
0.05

q

q

q
0.00

q       q     q                                            q

0.0          0.2       0.4            0.6        0.8       1.0

Values for theta

Here we can see that xobs = 8 has eﬀectively ruled out the values 0, 0.1, 0.2
and 1.0 for θ (0 and 1 have probability zero, the other values have very small
probabilities).

Continuous parameter

Now we consider the full range of possible values for θ, by taking Ω = [0, 1]. Because
Ω is continuous, the prior and the posterior are represented as probability density
functions. The main diﬀerence is that instead of summing to ﬁnd m, we have to
integrate:
m=           L(θ) p(θ) dθ.
Ω

We have to specify both the prior and the likelihood as functions of θ. Suppose,
initially, that the prior is uniform, i.e. that

1 θ ∈ [0, 1]
p(θ) =
0 otherwise.

In this case we have

3
>   punif <- function(theta) {
+       ifelse(0 <= theta, ifelse(theta <= 1, 1, 0))
+   }
>   L <- function(theta) {
+       stopifnot(all(0 <= theta & theta <= 1))
+       theta^8 * (1 - theta)^4
+   }

Both of these functions are written to work with a vector theta, which is useful for
drawing ﬁgures.
Now we can use the integrate function to ﬁnd m, and we can construct the
posterior probability density function:

> ii <- integrate(f = function(theta) {
+     L(theta) * punif(theta)
+ }, lower = 0, upper = 1)
> ii

0.0001554002 with absolute error < 1.7e-18

> m <- ii\$value
> post <- function(theta) {
+     L(theta) * punif(theta)/m
+ }

Note how we can write the function L(θ) × p(θ) directly into the f argument of
integrate.
Now we plot the result:

>   tt <- seq(from = 0, to = 1, length = 101)
>   plot(tt, post(tt), type = "l", lty = 1, lwd = 2,
+       main = "Prior and posterior probability density functions",
+       xlab = "Values of theta", ylab = "Probability density",
+       bty = "n")
>   lines(tt, punif(tt), lty = 2, lwd = 2)
>   legend("topleft", legend = c("Prior", "Posterior"),
+       lty = c(2, 1), lwd = 2, pch = NA, bty = "n")

4
Prior and posterior probability density functions

Prior
3.0

Posterior
2.5
Probability density

2.0
1.5
1.0
0.5
0.0

0.0         0.2     0.4        0.6          0.8         1.0

Values of theta

Q2. Repeat the analysis for a triangular prior:

> ptri <- function(theta) {
+     ifelse(0 <= theta, ifelse(theta <= 1, 2 *
+         (1 - theta), 0))
+ }

and contrast your result with that of the uniform prior.

Using the posterior

All inference about θ is based on the posterior distribution. For a point estimate, we
might choose the expectation or, less commonly, the median or the mode. Here are
examples of these calculations based on the uniform prior.
For the expectation, we need to ﬁnd

E{θ | X = xobs } =       θ q(θ) dθ,
Ω

so we will need to use integrate again.

5
> postExp <- integrate(f = function(theta) {
+     theta * post(theta)
+ }, lower = 0, upper = 1)\$value
> postExp

[1] 0.6428571

For the mode, we need our old friend optim:

> postmode <- optim(par = 0.5, fn = post, method = "L-BFGS-B",
+     lower = 0, upper = 1, control = list(fnscale = -1))\$par
> postmode

[1] 0.6666662

The median is slightly harder: it is that value θ† for which

θ†
q(θ) dθ = 0.5.
0

We can ﬁnd this by trial and error using integrate, but the more adventurous might
want to use uniroot:

>   uu <- uniroot(f = function(tdag) {
+       integrate(post, lower = 0, upper = tdag)\$value -
+           0.5
+   }, lower = 0, upper = 1)
>   postmedian <- uu\$root
>   postmedian

[1] 0.649848

Q3. Find the posterior standard deviation using integrate.

Q4. Find the posterior symmetric 95% credible interval using uniroot to compute
the 2.5th and 97.5th percentiles.

Q5. (hard) Find the 95% High Density Region (HDR)—this will require both uniroot
and integrate. Plot the result.

6

```
DOCUMENT INFO
Shared By:
Categories:
Stats:
 views: 6 posted: 2/6/2010 language: English pages: 6
How are you planning on using Docstoc?