Exercise 4.1

For the heavy-sleeper problem (see the Case study in Chapter 3), use the following priors to get posterior distribution (probabilities) in a future sample of \(n^*=20\) students.

  1. discrete prior
library(LearnBayes)
p <- seq(0.05, 0.95, by = 0.1)
prior <- c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.7, 0.1, 0, 0)
prior <- prior / sum(prior)
ns <- 20 ; ys <- 0:20
data <- c(11, 16)
post <- pdisc(p, prior, data)
pred <- pdiscp(post, prior, ns, ys)
cbind(round(pred, 3))
       [,1]
 [1,] 0.274
 [2,] 0.088
 [3,] 0.091
 [4,] 0.077
 [5,] 0.059
 [6,] 0.048
 [7,] 0.049
 [8,] 0.055
 [9,] 0.059
[10,] 0.059
[11,] 0.052
[12,] 0.040
[13,] 0.026
[14,] 0.014
[15,] 0.006
[16,] 0.002
[17,] 0.001
[18,] 0.000
[19,] 0.000
[20,] 0.000
[21,] 0.000
  1. Uniform prior
library(LearnBayes)
library(TailRank)
ab <- c(12, 17)
u <- 12 ; v <- 17
ns <- 20 ; ys <- 0:20
pprobs1 <- pbetap(ab, ns, ys)
pprobs2 <- dbb(ys, ns, u, v)
round(cbind(ys, pprobs1, pprobs2), 3)
      ys pprobs1 pprobs2
 [1,]  0   0.000   0.000
 [2,]  1   0.003   0.003
 [3,]  2   0.010   0.010
 [4,]  3   0.025   0.025
 [5,]  4   0.049   0.049
 [6,]  5   0.078   0.078
 [7,]  6   0.108   0.108
 [8,]  7   0.129   0.129
 [9,]  8   0.137   0.137
[10,]  9   0.131   0.131
[11,] 10   0.112   0.112
[12,] 11   0.086   0.086
[13,] 12   0.059   0.059
[14,] 13   0.037   0.037
[15,] 14   0.020   0.020
[16,] 15   0.009   0.009
[17,] 16   0.004   0.004
[18,] 17   0.001   0.001
[19,] 18   0.000   0.000
[20,] 19   0.000   0.000
[21,] 20   0.000   0.000
  1. beta prior
a <- 3.26 ; b <- 7.19
s <- 11 ; f <- 15
ab <- c(a + s, b + f)
ns <- 20 ; ys <- 0:20
pred1 <- pbetap(ab, ns, ys)
pred2 <- dbb(ys, ns, a + s, b + f)
round(cbind(ys, pred1, pred2), 3)
      ys pred1 pred2
 [1,]  0 0.000 0.000
 [2,]  1 0.003 0.003
 [3,]  2 0.012 0.012
 [4,]  3 0.031 0.031
 [5,]  4 0.059 0.059
 [6,]  5 0.093 0.093
 [7,]  6 0.123 0.123
 [8,]  7 0.142 0.142
 [9,]  8 0.144 0.144
[10,]  9 0.128 0.128
[11,] 10 0.102 0.102
[12,] 11 0.072 0.072
[13,] 12 0.045 0.045
[14,] 13 0.025 0.025
[15,] 14 0.012 0.012
[16,] 15 0.005 0.005
[17,] 16 0.002 0.002
[18,] 17 0.000 0.000
[19,] 18 0.000 0.000
[20,] 19 0.000 0.000
[21,] 20 0.000 0.000

Exercise 4.2

Use the simulation method to get the predictive distribution above.

  1. discrete prior
m<-2000
p<-seq(0.05,0.95,by = 0.1)
prior <- c(1, 5.2, 8, 7.2, 4.6, 2.1, 0.7, 0.1, 0, 0)
prior <-prior/sum(prior)
data<-c(11,16)
post<-pdisc(p,prior,data)
p_samples<-sample(p,size = m, replace = TRUE,prob = post)
#sample takes a sample of the specified size from the elements of x .
ns<-20
y_samples<-rbinom(m,ns,p_samples)
freq<-table(y_samples)
ys<-as.integer(names(freq))
predprob<-freq/sum(freq)
round(cbind(predprob),3)
   predprob
0     0.001
1     0.004
2     0.014
3     0.035
4     0.063
5     0.106
6     0.119
7     0.142
8     0.151
9     0.126
10    0.095
11    0.068
12    0.040
13    0.020
14    0.010
15    0.004
16    0.002
17    0.000
plot(ys,predprob,type = "h",xlab = "y",
     ylab = "Predictive Probability with 
     discrete prior",cex.lab =1, cex.axis =1.0)
pred<-pdiscp(p,post,ns,ys)
points(ys,pred,type = "p", col = "red")

  1. Uniform prior
m = 1000 
ab <- c(12, 17)
u <- 12 ; v <- 17
ns <- 20 ; ys <- 0:20
p = rbeta(m, u, v) 
y = rbinom(m, ns, p) 
freq = table(y) 
ys = as.integer(names(freq)) 
predprob = freq/sum(freq) 
plot(ys, predprob, type="h", xlab="y", 
     ylab="Predictive Probability with Uniform prior", 
     cex.lab=1, cex.axis=1.0) 
upred1 = pbetap(ab, ns, ys) 
upred2 = dbb(ys, ns, u, v) 
points(ys, upred1, type="p", pch='o', col='blue')
points(ys, upred2, type="p", pch='*', col='red')

  1. Beta prior
m = 1000 
a = 3.26; b = 7.19 
s = 11; f = 16 
ns = 20 
p = rbeta(m, a+s, b+f) 
y = rbinom(m, ns, p) 
freq = table(y) 
ys = as.integer(names(freq)) 
predprob = freq/sum(freq) 
plot(ys, predprob, type="h", xlab="y", 
     ylab="Predictive Probability", 
     cex.lab=1, cex.axis=1.0) 
library(LearnBayes) 
TruePred1 = pbetap(c(a+s,b+f), ns, ys) 
points(ys, TruePred1, type="p", col='red') 
library(TailRank) 
TruePred2 = dbb(ys, ns, a+s, b+f) 
points(ys, TruePred2, type="p", pch='*', col='blue')

Exercise 5.1 — Problem 5.3, Cowles (2013)

Verify analytically that a uniform prior on π induces the density in (5.2) as the prior on \(\phi = logit(\pi)\)

\[\begin{aligned} \pi \sim U(0,1)\Rightarrow P_{\pi}(p)=1,0<p<1\\ \phi=g(p)=logit(p)=log\frac{p}{1-p}\\ p=g^{-1}(\phi)·\frac{e^{\phi}}{1+e^{\phi}}=:h(\phi)\\ P_{\phi}(\phi)=P_{\pi}(h{\phi})·h'(\phi)=1·(\frac{e^{\phi}}{1+e^{\phi}})'=\frac{e^{\phi}}{(1+e^{\phi})^2} \end{aligned}\]

Exercise 5.2 — Problem 5.4, Cowles (2013)

Suppose that a trucking company owns a large fleet of well-maintained trucks and assume that breakdowns appear to occur at random times. The president of the company is interested in learning about the daily rate \(\lambda\) at which breakdowns occur. (Realistically, each truck would have a breakdown rate that depends possibly on its type, age, condition, driver, usage, etc. The breakdown rate for the whole company can be viewed as the sum of the breakdown rates of the individual trucks.) For a given value of the rate parameter \(\lambda\) , it is known that the number of breakdowns y on a particular day has a Poisson distribution with mean \(\lambda\):

\[p(y|\lambda)=\frac{e^{-\lambda}\lambda^y}{y!}\]

  1. Suppose that one observes the number of truck breakdowns for n consecutive days, denoted as \(y_1, . . . , y_n\). If one assumes that these are exchangeable measurements (conditionally independent given \(\lambda\)), find the joint probability distribution of \(y_1, . . . , y_n\)

\(Sol.\)

\[P(Y_1=y_1,Y_2=y_2,...,Y_n=y_n|\lambda)=P(y_1|\lambda)···P(y_n|\lambda)=e^{-\lambda n}\frac{\lambda^{\sum_{i=1}^{n}y_i}}{\Pi_{i=1}^{n}y_i!}\]

  1. The numbers of breakdowns for 5 days are recorded to be 2, 5, 1, 0, and 3. Find the likelihood function L(\(\lambda\)) of the rate parameter \(\lambda\) for these observations. Graph this function. (You may either use R or do it “by hand” by calculating the likelihood for the values from 0 to 10 by step 0.1 and connecting the points with a smooth curve.)

\(Sol.\)

\[L(\lambda)=e^{-\lambda n}\frac{\lambda^{\sum_{i=1}^{n}y_i}}{\Pi_{i=1}^{n}y_i!}\]

lambda <- seq(0, 10, 0.1)
n <- 5
y <- c(2, 5, 1, 0, 3)
y_sum <- sum(y)
y_fac <- prod(factorial(y))
likelihood <- exp(-lambda * n) * lambda^y_sum / y_fac
log_likelihood <- log(likelihood)
par(mfrow = c(2,1), mar = c(4,2,1,1))
plot(lambda, likelihood, type = "l", main = "likelihood")
plot(lambda, log_likelihood, type = "l", main = "log likelihood")

  1. Use calculus to find the MLE of \(\lambda\). Then use the \(poisson.test\) function in R to confirm the MLE and to obtain a 95% frequentist confidence interval for \(\lambda\).

\[\begin{aligned} &L(\lambda)=e^{-\lambda n}\frac{\lambda^{\sum_{i=1}^{n}y_i}}{\Pi_{i=1}^{n}y_i!}\\ &\frac{\partial L(\lambda)}{\partial \lambda}=\frac{1}{\Pi_{i=1}^{n}y_i!}[-ne^{-\lambda n}\lambda^{\sum_{i=1}^{n}y_i}+e^{-\lambda n}\sum_{i=1}^{n}y_i\lambda^{(\sum_{i=1}^{n}y_i)-1}]=0\\ &-n+\frac{\sum_{i=1}^{n}y_i}{\lambda}=0,n=5,\sum_{i=1}^{n}y_i=11\\ &\lambda = 2.2 \end{aligned}\]

Thus the MLE is 2.2.

# Use poisson.test function to obtain MLE
poisson_test_result <- poisson.test(x = y_sum, T = length(y),r=2.2)
# Print the result
print(poisson_test_result)

    Exact Poisson test

data:  y_sum time base: length(y)
number of events = 11, time base = 5, p-value = 1
alternative hypothesis: true event rate is not equal to 2.2
95 percent confidence interval:
 1.098232 3.936408
sample estimates:
event rate 
       2.2

Then we have the MLE \(\hat{\lambda}=2.2\) and a 95% frequentist confidence interval for \(\lambda\) [1.098232, 3.936408].

Exercise 5.3 — Problem 5.5, Cowles (2013)

The president of the company has some knowledge about the location of the Poisson rate parameter \(\lambda\) based on the observed number of breakdowns from previous years. His prior beliefs about \(\lambda\) are represented in the following:

\[p(\lambda)\propto\lambda^3exp(-2\lambda),\lambda>0\]

  1. Is this prior a member of a particular parametric family? If so, what family and what are the prior parameters?

    From Gamma distribution famiy. \(\alpha = 4,\beta= 2\)

  2. Plot this prior density in R. Based on the plot, describe the president’s prior beliefs about the rate parameter \(\lambda\).

# Define the prior density function
prior_density <- function(lambda) {
  return(lambda^3 * exp(-2*lambda))
}
# Generate a sequence of lambda values
lambda_values <- seq(0, 10, by = 0.01)
# Calculate the prior density values
density_values <- prior_density(lambda_values)
# Calculate the integral of the density function
integral <- integrate(prior_density, lower = 0, upper = Inf)$value
# Normalize the density function
normalized_density <- density_values / integral
# Plot the normalized prior density
plot(lambda_values, normalized_density, type = "l", 
     main = "Prior Density for Rate Parameter λ",
     xlab = "λ", ylab = "Density")

  1. Write out the mathematical form of the unnormalized posterior density. Identify its parametric family and parameters.

\[\begin{aligned} p(\lambda \mid x) &\propto p(x \mid \lambda) p(\lambda)\\ &=e^{-\lambda n}\frac{\lambda^{\sum_{i=1}^{n}y_i}}{\Pi_{i=1}^{n}y_i!}\lambda^3e^{-2\lambda}\\ &=e^{-(n+2)\lambda}\lambda^{3+\sum_{i=1}^{n}y_i}\frac{1}{\Pi_{i=1}^{n}y_i!}\\ &=\frac{1}{1440}\lambda^{14}e^{-7\lambda} \end{aligned}\]

\Gamma(\alpha,\lambda)=\Gamma(15,7)

  1. Find the posterior mean and 95% central credible set for \(\lambda\) based on this posterior.

    The posterior mean is \(\frac{15}{7}=2.142857\)

# Prior parameters
alpha <- 4
lambda <- 2
# Likelihood parameters
n <- length(y)
sum_y <- sum(y)
# Posterior parameters
alpha_posterior <- alpha + sum_y
lambda_posterior <- lambda + n
# Posterior mean
posterior_mean <- alpha_posterior / lambda_posterior
# Posterior 95% credible interval
lower <- qgamma(0.025, shape = alpha_posterior, rate = lambda_posterior)
upper <- qgamma(0.975, shape = alpha_posterior, rate = lambda_posterior)
# Print results
cat("Posterior mean of λ:", posterior_mean, "\n")
Posterior mean of λ: 2.142857 
cat("95% central credible set for λ: (", lower, ",", upper, ")\n")
95% central credible set for λ: ( 1.199341 , 3.35566 )

  1. Was the president’s prior from a conjugate family for the Poisson likelihood? How could you tell?

    Yes, the president’s prior for the Poisson likelihood was from a conjugate family. We can tell because the prior distribution (gamma distribution) belongs to a conjugate family for the Poisson likelihood.

Exercise 6.2 — Problem 6.2, Cowles (2013)

The observed weights (in grams) of 20 pieces of candy randomly sampled from candy-making machines in a certain production area are as follows:

46, 58, 40, 47, 47, 53, 43, 48, 50, 55,

49, 50, 52, 56, 49, 54, 51, 50, 52, 50.

Assume that weights of this type of candy are known to follow a normal distribution, and that the mean weight of candies produced by machines in this area is known to be 51 g. We are trying to estimate the variance, which we will now call θ.

  1. What is the conjugate family of prior distributions for a normal variance (not precision) when the mean is known?

Conjugate prior for a normal variance when tne mean is known: inverse gamma \(IG(\alpha, \beta)\)

  1. Suppose previous experience suggests that the expected value of θ is 12 and the variance of θ is 4. What parameter values are needed for the prior distribution to match these moments? We should get \(\alpha,\beta\) through trial and error.

\(\alpha=38,\beta=444\)

  1. What is the posterior distribution π(θ|y) for these data under the prior from the previous step? \[\begin{aligned} \pi(\theta|y) \propto \frac{1}{\theta^{n/2+\alpha+1}}exp[-\frac{1}{\theta}(\frac{nv}{1}+\beta)],0<\theta<\infty\\ i.e. \theta|y \sim IG(\alpha + \frac{n}{2}, \beta + \frac{nv}{2})=IG(48, 628) \end{aligned}\]
alpha = 38
beta = 444
mu = 51
y = c(46,58,40,47,47,53,43,48,50,55,49,50,52,56,49,54,51,50,52,50)
n = 20
v = sum((y-mu)^2)/n
#posterior parameter
alpha1 = alpha + n/2 
beta1 = beta + v*n/2 
c(alpha1, beta1)
[1]  48 628
  1. Find the posterior mean and variance of θ.

\[ E(\theta \mid y)=\frac{\beta}{\alpha-1}=\frac{628}{47}=13.3617\\ Var(\theta \mid y)=\frac{\beta^2}{(\alpha-1)^2(\alpha-2)}=\frac{628^2}{47^246}=3.881197\\ \]

Exercise 6.2 — Exercise 3.2, Albert (2009)

Suppose a random sample \(y_1,y_2,...,y_n\) is taken from an exponential distribution with mean \(theta\)

\(p(y|\theta)=\theta^{-1}exp\{-y/\theta\}\)

  1. Show that if we assign the usual noninformative prior \(\pi(\theta) \propto 1/\theta\), then the posterior density is given, up to a proportionality constant, by

\(\pi(\theta|data) \propto \theta^{-n-1}exp\{-s/\theta\}\),

where \(s=\sum_{i=1}^{n}y_i\)

\(Proof:\)

The likelihood function is:

\(L(\theta|y)=\Pi_{i=1}^{n}\theta^{-1}exp\{-y/\theta\}=\theta^{-n}exp\{-\sum_{i=1}^{n}y_i/\theta\}\);

The prior is:

\(\pi(\theta) \propto 1/\theta\);

Then the posterior density can be written as:

\(\pi(\theta|data) \propto L(\theta|y) \times \pi(\theta) \propto \theta^{-n-1}exp\{-\sum_{i=1}^{n}y_i/\theta\}=\theta^{-n-1}exp\{-s/\theta\}\).

  1. Show that if we transform \(\theta\) to \(\lambda = 1/\theta\), then \(\lambda\) has a gamma density with shape parameter \(n\) and rate parameter \(s\).

\(Proof:\)

If we transform \(\theta\) to \(\lambda = 1/\theta\), \(\pi(\lambda|data) \propto \lambda^{n+1}exp\{-s\lambda\}\sim Ga(n,s)\)

  1. In a life-testing illustration, five bulbs are tested with observed burn times (in hours) of 751, 594, 1213, 1126 and 819. Using the R function rgamma(), simulate 1000 values from the posterior distribution of \(\lambda\).
burntime <- c(751, 594, 1213, 1126, 819)
s <- sum(burntime)
n <- 5
lambda <- rgamma(1000, shape = n, rate = s)
plot(density(lambda), main = "Buring times", xlab = "lambda")

  1. By transforming these simulated draws, obtain a simulated sample from the posterior distribution of \(\theta\).
theta_samples <- 1 / lambda
plot(density(theta_samples), main = "Posterior distribution of theta", xlab = "theta")

  1. Estimate the posterior probability that the mean life time \(\theta\) exceeds 1000 hours.
sum(theta_samples > 1000) / length(theta_samples)
[1] 0.468