Exercise 4.1

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

Solution

  1. discret 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)
data<-c(11,16)#sample information
post<-pdisc(p,prior,data)
ns = 20
ys = 0:20
pred = pdiscp(p,post,ns,ys)
round(cbind(ys,pred),3)
      ys  pred
 [1,]  0 0.001
 [2,]  1 0.004
 [3,]  2 0.014
 [4,]  3 0.035
 [5,]  4 0.065
 [6,]  5 0.099
 [7,]  6 0.129
 [8,]  7 0.145
 [9,]  8 0.143
[10,]  9 0.125
[11,] 10 0.096
[12,] 11 0.066
[13,] 12 0.040
[14,] 13 0.022
[15,] 14 0.010
[16,] 15 0.004
[17,] 16 0.001
[18,] 17 0.000
[19,] 18 0.000
[20,] 19 0.000
[21,] 20 0.000
#max(pred)
which.max(pred)
[1] 8
#which(abs(pred-0.02030242)<1e-5)# 寻找任何一个数据所在位置
########################(2)uniform prior
library(LearnBayes)
a=1;b=1#U(0,1)=Beta(1,1)
s=11; f=16#sample information
ab=c(a+s,b+f)
ns = 20
ys = 0:20
pred2_1 = pbetap(ab, ns, ys) 
round(cbind(ys,pred2_1), 3)
      ys pred2_1
 [1,]  0   0.000
 [2,]  1   0.003
 [3,]  2   0.010
 [4,]  3   0.025
 [5,]  4   0.049
 [6,]  5   0.078
 [7,]  6   0.108
 [8,]  7   0.129
 [9,]  8   0.137
[10,]  9   0.131
[11,] 10   0.112
[12,] 11   0.086
[13,] 12   0.059
[14,] 13   0.037
[15,] 14   0.020
[16,] 15   0.009
[17,] 16   0.004
[18,] 17   0.001
[19,] 18   0.000
[20,] 19   0.000
[21,] 20   0.000
#max(pred2)
which.max(pred2_1)
[1] 9
#方法2
library(TailRank)
a=1;b=1#U(0,1)=Beta(1,1)
s=11; f=16#sample information
ns = 20
ys = 0:20
pred2_2 <-dbb(ys,ns,a+s,b+f) 
round(cbind(ys,pred2_1,pred2_2), 3)
      ys pred2_1 pred2_2
 [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
########################(3)beta prior
library(LearnBayes)
a=3.26;b=7.19
s=11; f=16
ab=c(a+s,b+f)
ns = 20#new sample size n*=20
ys = 0:20#y* successes 
pred3 = pbetap(ab, ns, ys) 
round(cbind(ys,pred3), 3)
      ys pred3
 [1,]  0 0.001
 [2,]  1 0.004
 [3,]  2 0.014
 [4,]  3 0.035
 [5,]  4 0.066
 [6,]  5 0.100
 [7,]  6 0.130
 [8,]  7 0.145
 [9,]  8 0.143
[10,]  9 0.124
[11,] 10 0.096
[12,] 11 0.065
[13,] 12 0.040
[14,] 13 0.021
[15,] 14 0.010
[16,] 15 0.004
[17,] 16 0.001
[18,] 17 0.000
[19,] 18 0.000
[20,] 19 0.000
[21,] 20 0.000
#max(pred3)
which.max(pred3)
[1] 8
#方法2
a<-3.26;b<-7.19
s<-11;f<-16
ns<-20;ys<-0:20
pred3_1<-dbb(ys,ns,a+s,b+f)
round(cbind(ys,pred3,pred3_1), 3)
      ys pred3 pred3_1
 [1,]  0 0.001   0.001
 [2,]  1 0.004   0.004
 [3,]  2 0.014   0.014
 [4,]  3 0.035   0.035
 [5,]  4 0.066   0.066
 [6,]  5 0.100   0.100
 [7,]  6 0.130   0.130
 [8,]  7 0.145   0.145
 [9,]  8 0.143   0.143
[10,]  9 0.124   0.124
[11,] 10 0.096   0.096
[12,] 11 0.065   0.065
[13,] 12 0.040   0.040
[14,] 13 0.021   0.021
[15,] 14 0.010   0.010
[16,] 15 0.004   0.004
[17,] 16 0.001   0.001
[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.

不对：

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)
p_samples<-sample(p,size = m, replace = TRUE,prob = prior)
#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.020
1     0.050
2     0.077
3     0.090
4     0.095
5     0.115
6     0.108
7     0.099
8     0.089
9     0.070
10    0.054
11    0.044
12    0.032
13    0.030
14    0.013
15    0.007
16    0.003
17    0.003
18    0.001
plot(ys,predprob,type = "h",xlab = "y",
     ylab = "Predictive Probability with 
     discrete prior",cex.lab =1, cex.axis =1.0)
pred<-pdiscp(post,prior,ns,ys)
points(ys,pred,type = "p", col = "red")

正确

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
1     0.004
2     0.018
3     0.039
4     0.064
5     0.105
6     0.121
7     0.153
8     0.138
9     0.126
10    0.097
11    0.054
12    0.040
13    0.025
14    0.009
15    0.003
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")

下面这段代码不对

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,post)
freq<-table(y_samples)
ys<-as.integer(names(freq))
predprob<-freq/sum(freq)
round(cbind(predprob),3)
   predprob
0     0.634
1     0.066
2     0.044
3     0.040
4     0.022
5     0.032
6     0.027
7     0.026
8     0.026
9     0.025
10    0.020
11    0.016
12    0.011
13    0.005
14    0.002
15    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")

- 2.Uniform prior

m = 1000 
a= 1; b = 1 
s = 11; f = 16 
ns = 20 
p = rbeta(m, a+s, b+f) #生成后验分布横坐标p的随机数1000
y = rbinom(m, ns, p) # 在ns次试验中随机生成m个数，其中每一个数表示成功的次数。
freq = table(y) # 列表y中的数字出现的次数
ys = as.integer(names(freq)) 
predprob = freq/sum(freq) 
round(cbind(predprob),3)
   predprob
0     0.001
1     0.004
2     0.009
3     0.025
4     0.056
5     0.060
6     0.101
7     0.140
8     0.121
9     0.139
10    0.114
11    0.097
12    0.061
13    0.037
14    0.019
15    0.013
16    0.003
plot(ys, predprob, type="h", xlab="y", ylab="Predictive Probability", cex.lab=2, cex.axis=2.0)

TruePred2 = pbetap(c(a+s,b+f), ns, ys)
points(ys, TruePred2, type="p", col='red')

m = 1000 
a3 = 3.26; b3 = 7.19 
s = 11; f = 16 
ns = 20 
p = rbeta(m, a3+s, b3+f) #生成后验分布横坐标p的随机数1000
y = rbinom(m, ns, p) # 在ns次试验中随机生成m个数，其中每一个数表示成功的次数。
freq = table(y) # 列表y中的数字出现的次数
ys = as.integer(names(freq)) 
predprob = freq/sum(freq) 
round(cbind(predprob),3)
   predprob
1     0.003
2     0.015
3     0.028
4     0.055
5     0.093
6     0.117
7     0.170
8     0.156
9     0.133
10    0.094
11    0.065
12    0.034
13    0.022
14    0.012
15    0.003
plot(ys, predprob, type="h", xlab="y", ylab="Predictive Probability", cex.lab=2, cex.axis=2.0)

TruePred3 = pbetap(c(a3+s,b3+f), ns, ys)
points(ys, TruePred3, type="p", col='red')

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 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\mid\lambda)=\frac{e^{-\lambda}\lambda^y}{y!},y=0,1,2,\cdots \]

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

Since the breakdowns are assumed to be independent and identically distributed (i.i.d) with a Poisson distribution parameterized by \(\lambda\), the joint probability distribution of \(y_1, \dots, y_{n}\) is the product of the individual probability distributions: \[ p\left(y_1, \ldots, y_n \mid \lambda\right)=\prod_{i=1}^n p\left(y_i \mid \lambda\right)=\prod_{i=1}^n \frac{e^{-\lambda} \lambda^{y_i}}{y_{i} !}=e^{-\lambda n}\frac{\lambda^{\sum_{i=1}^{n}y_i}}{\Pi_{i=1}^{n}y_i!} \]

  1. The number 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)

\[ 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")

The likelihood function \(\lambda\) is the probability of observing the given breakdown counts for 5 days, given the parameter \(\lambda)\). Given the breakdown counts $ 2,5,1,0 $, and \(3\), the likelihood function is: \[ L(\lambda)=p(2,5,1,0,3 \mid \lambda)=\frac{e^{-5 \lambda} \lambda^{11}}{2 ! \cdot 5 ! \cdot 1 ! \cdot 0 ! \cdot 3 !} \]

# define the likelihood function
likelihood <- function(lambda){
  numerator <- exp(-5 * lambda) * lambda^11
  denominator <- factorial(2) * factorial(5) * factorial(1) * factorial(0) * factorial(3)
  return(numerator / denominator)
}

# calculate the likelihood for values from 0 to 10 by step 0.1
lambda_seq <- seq(0, 10, by = 0.1)
likelihood_seq <- sapply(lambda_seq, likelihood)
#likelihood_seq 最终会是一个向量，其中包含了对 lambda_seq 序列中每个值应用 likelihood 函数的结果
# plot the likelihood curve
plot(lambda_seq, likelihood_seq, type = "l", xlab = "lambda", ylab = "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\).

To find the maximum likelihood estimate (MLE) of \(\lambda)\), we need to maximize the likelihood function \(L(\lambda)\) with respect to \(\lambda\). Taking the natural logarithm of \(L(\lambda)\) (which doesn’t change the location of the maximum), we have: \[ \ln L(\lambda)=-5 \lambda+11 \ln (\lambda)+C \] where \(C\) is a constant term that doesn’t affect the location of the maximum. To find the MLE, we take the derivative of \(\ln L(\lambda)\) with respect to \(\lambda\), set it equal to zero, and solve for \(\lambda\) : \[ \begin{aligned} & \frac{d}{d \lambda} \ln L(\lambda)=-5+\frac{11}{\lambda}=0 \\ & \Rightarrow \frac{11}{\lambda}=5 \\ & \Rightarrow \lambda=\frac{11}{5} \end{aligned} \] - 方法2

\[ \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} \]

# 定义数据
data <- c(2, 5, 1, 0, 3)

# 计算总和
sum_data <- sum(data)

# 计算MLE
lambda_mle <- sum_data / length(data)
lambda_mle
[1] 2.2
# 使用poisson.test函数进行推断
poisson_test_result <- poisson.test(x = sum_data, T = 5, 2.2,alternative = "two.sided", conf.level = 0.95)
print(poisson_test_result)

    Exact Poisson test

data:  sum_data time base: 5
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^3 \exp (-2 \lambda), \quad \lambda>0 \]

  1. Is this prior a member of a particular parametric family?If so, what family and what are the prior parameters?
  1. Plot this prior density in R. Based on the plot,describe the president’s prior beliefs about the rate parameter \(\lambda\).
par(mfrow=c(2,1),mar=c(4,2,1,1))
#unnormalized
lambda <- seq(0, 10, by=0.01)
prior_density <- lambda^3 * exp(-2*lambda)
plot(lambda, prior_density, type='l', xlab='lambda', ylab='Density', main='Prior Density Plot')

#normalized
prior_density<-function(lambda){
  return(lambda^3 * exp(-2*lambda))
}
  lambda <- seq(0, 10, by=0.01)
  density_values<-prior_density(lambda)
  integral<-integrate(prior_density,lower =0,
                      upper = Inf)$value
  normalized_density<-density_values/integral
  plot(lambda,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.

后验密度的数学形式为: \[ p(\lambda \mid x) \propto p(x \mid \lambda) p(\lambda) \propto \frac{e^{-5 \lambda} \lambda^{11}}{2 ! \cdot 5 ! \cdot 1 ! \cdot 0 ! \cdot 3 !} \times \lambda^3 \exp (-2 \lambda) \] - 这个后验分布的形式还是Gamma分布，但参数会有所变化。\(\alpha=15， \beta=7\) 。 \[\begin{aligned} \pi(\lambda|X) &\propto L(\lambda) \times prior\\ &=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}\]

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

\[ E(\lambda \mid x)=\frac{\alpha}{\beta}=\frac{15}{7}\\ Var(\lambda \mid x)=\frac{\alpha}{\beta^2}=\frac{15}{7^2} \]

#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
qgamma(c(0.025,0.975),alpha_posterior,lambda_posterior)
[1] 1.199341 3.355660
#print results
cat("posterior mean of lambda:",posterior_mean)
posterior mean of lambda: 2.142857
  1. Was the president’s prior from a conjugate family for the Poisson likelihood? How could you tell?

Exercise 6.1——-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:

weights<-c(46, 58, 40, 47, 47, 53, 43, 48, 50 ,55, 49, 50, 52, 56, 49, 54, 51, 50, 52, 50)
mean(weights)
[1] 50
cbind(weights)
      weights
 [1,]      46
 [2,]      58
 [3,]      40
 [4,]      47
 [5,]      47
 [6,]      53
 [7,]      43
 [8,]      48
 [9,]      50
[10,]      55
[11,]      49
[12,]      50
[13,]      52
[14,]      56
[15,]      49
[16,]      54
[17,]      51
[18,]      50
[19,]      52
[20,]      50
x_i<-(weights-52)^2
mean(x_i)
[1] 21.4

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 \mathrm{~g}\). We are trying to estimate the variance, which we will now call \(\theta\).

-solution

\[ p(y_1,y_2,\cdots,y_n\mid\mu,\sigma^2)=L(y\mid\sigma^2)\\=\Pi_{i=1}^n(2\pi\sigma^2)^{-\frac{1}{2}}\exp\{{-\frac{\sum_{i=1}^n(y_i-\mu)^2}{2\sigma^2}}\}\\ =(2\pi\sigma^2)^{-\frac{n}{2}}\exp\{{-\frac{nv}{2\sigma^2}}\}\\,where\quad v=\frac{1}{n}\sum_{i=1}^n(y_i-\mu)^2 \] \[ L(y\mid\sigma^2)\propto(1/\sigma^2)^{\frac{n}{2}}\exp\{{-\frac{nv}{2\sigma^2}}\}\\ y\mid\sigma^2\sim IG(\frac{n}{2}-1,\frac{nv}{2}) \]

\[ \sigma^2\sim IG(\alpha,\beta)\\ \pi(\sigma^2)=\frac{\beta^\alpha}{\Gamma(\alpha)} {(\frac{1}{\sigma^2})}^{\alpha+1}\exp(-\frac{\beta}{\sigma^2})\\ \sigma^2\sim IG(\frac{n_0}{2},\frac{n_0\sigma_{0}^2}{2}) \]

\[ p(\sigma^2\mid y)=L(y\mid\sigma^2)\pi(\sigma^2)\\ =(2\pi\sigma^2)^{-\frac{n}{2}}\exp\{{-\frac{nv}{2\sigma^2}}\}*\frac{\beta^\alpha}{\Gamma(\alpha)} {(\frac{1}{\sigma^2})}^{\alpha+1}\exp(-\frac{\beta}{\sigma^2})\\ \propto {(\frac{1}{\sigma^2})}^{\frac{n}{2}+\alpha+1}\exp(-\frac{\beta+\frac{nv}{2}}{\sigma^2}) \] 即 \[ \sigma^2\mid y\sim IG(\frac{n}{2}+\alpha,\beta+\frac{nv}{2}) \] - (b)Suppose previous experience suggests that the expected value of \(\theta\) is 12 and the variance of \(\theta\) is 4 . What parameter values are needed for the prior distribution to match these moments?

\[ E(\theta)=\frac{\beta}{\alpha-1}=12\\ Var(\theta)=\frac{\beta^2}{(\alpha-1)^2(\alpha-2)}=4\\ \alpha=38,\beta=444\\ \theta\sim IG(38,444) \]

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

\[ \theta \mid y\sim IG(48,628) \] - d. Find the posterior mean and variance of \(\theta\).

\[ 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, \ldots, y_n\) is taken from an exponential distribution with mean \(\theta\) \[ p(y \mid \theta)=\theta^{-1} \exp \{-y / \theta\} . \] (a) 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 \mid \text { data }) \propto \theta^{-n-1} \exp \{-s / \theta\} \]

\[ L(y \mid \theta)=p(y \mid \theta)=\\ \Pi_{i=1}^n\theta^{-1} \exp \{-y_i / \theta\} =\frac{1}{\theta^{n}}\exp \{-\sum_{i=1}^n y_i/ \theta\}\\=\frac{1}{\theta^{n}}\exp \{-s/ \theta\}. \]

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

\[ \pi(\theta \mid \text { data })=p(y \mid \theta)\pi(\theta)\propto\frac{1}{\theta^{n+1}}\exp \{-s/ \theta\}. \]

burntime<-c(751,594, 1213, 1126, 819)
sum<-sum(burntime)
n<-5
lambda<-rgamma(1000,shape=n,rate=sum)
plot(density(lambda),main= "burn times",xlab= "lambda")

  1. By transforming these simulated draws, obtain a simulated sample from the posterior distribution of \(\theta\).
rate<-c(751,594, 1213, 1126,819)
theta_samples<-1/rgamma(1000,5,sum(rate))
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.479
length(theta_samples)
[1] 1000
sample <- c(751, 594, 1213, 1126, 819)
n <- 5
s <- sum(sample)
# 计算theta <= 1000的概率，再用1减去
prob1 <- 1 - pinvgamma(1000, shape = n, rate = s)
# 也可以计算lambda <= 0.001的概率
prob2 <- pgamma(0.001, shape = n, rate = s)
print(cbind(prob1,prob2))