Given the hierarchical conjugate gamma-Poisson model described in the BUGS examples documentation, let’s derive a Gibbs Sample estimation to find out which pumps are more reliable:

The number of of failures \(x_i\) is assumed to follow a Poisson:

\[X_i \sim \textrm{Poisson}(\theta_i t_i) \hspace{5mm} i = 1,\dots,10\]

A conjugate gamma prior distribution is adopted for the failure rates:

\[\theta_i \sim \textrm{Gamma}(\alpha,\beta)\] And we also have the following prior specification for the hyperparameter \(\beta\): \[\alpha=1\] \[\beta \sim \textrm{Gamma}(0.1,1.0)\] The full and partial posterior conditionals are then:

\[h(\underline{\mathbf{\theta}}, \alpha, \beta| \underline{\mathbf{x}}, \underline{\mathbf{t}} ) = g(\underline{\mathbf{x}}, \underline{\mathbf{t}}| \underline{\mathbf{\theta}},\alpha, \beta)f_\theta(\underline{\mathbf{\theta}}|\alpha,\beta)f_\alpha(\alpha)f_\beta(\beta) \] Which then simplifies for each \(\theta_i\) and \(\beta\) to known distributions which can be Gibbs sampled:

\[\theta_j|\alpha,\beta,x_j,t_j) \sim \Gamma(x_j+1,t_j+\beta)\]

\[\beta|\underline{\mathbf{\theta}} \sim \Gamma(10\alpha+.1,\sum_{j=1}^{10}\theta_j+1)\]

Below is the Gibbs Sampling routine. Later, we shall also place a prior on \(\alpha\) which will turn into an unknown distribution for \(f_\alpha\) which we will have to implement through a M-H algirithm.

# Data
t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)
x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22)

n = 10; N = 10000; B = 1000
a = 1.0; b = 0.1; c = 1.0
alpha = beta  =numeric(N)
lower = upper = M = numeric(n)

# Initializing parameters
theta = matrix((x*t/sum(t)), nrow=N, ncol=n)
thet.a = matrix(0,nrow=N-B,ncol = n)
#alpha[1] =                       
beta = alpha = rep(1,N)

# The Gibbs Sampler, proper:
for(i in 2:N) {
  for(j in 1:n){   
    theta[i,j] = rgamma(1, shape = x[j] + alpha[1], rate = t[j] + beta[i-1])
  }
  beta[i] = rgamma(1, shape = b + n*alpha[1],  rate = c + sum(theta[i,]))
}

Let’s produce some standard plots for the parameters:

par(mfrow=c(3,4))
plot(beta[B:N],pch=".", main=expression(beta))
for(j in 1:3) plot(theta[B:N,j],pch=".",ylab=paste("theta[",j,"]"), main=expression(theta[j]))
hist(beta[B:N], probability=TRUE, main="")
for(j in 1:3) hist(theta[B:N,j],probability=TRUE,pch=".",main=NULL)
acf(beta[B:N], main="")
for(j in 1:3) acf(theta[B:N,j],main="")

credible intervals for \(\beta\) and \(\underline{\mathbf{\theta}}\)

library(forestplot)
## Loading required package: grid
## Loading required package: magrittr
## Loading required package: checkmate
# First the point estimates
nn = N-B; L = as.integer(.025*nn); U = as.integer(.975*nn)
bet.a=sort(beta[(B+1):N])
m.b = mean(bet.a);beta.l = bet.a[L]; beta.u = bet.a[U]
print(paste("Point & CI estimates for beta:",round(m.b,digits=4),"& (",round(beta.l,digits=4),",",round(beta.u,digits=4),")"))
## [1] "Point & CI estimates for beta: 1.3428 & ( 0.5757 , 2.495 )"
for (j in 1:n){
  thet.a[,j] = sort(theta[(B+1):N,j])
  M[j] = mean(thet.a[,j])
  lower[j] = thet.a[L,j]
  upper[j] = thet.a[U,j]
  print(paste("Point & CI estimates for theta(",j,"):",round(M[j],digits=4),"& (",round(lower[j],digits=4),",",round(upper[j],digits=4),")"))
  
}
## [1] "Point & CI estimates for theta( 1 ): 0.0627 & ( 0.0232 , 0.1223 )"
## [1] "Point & CI estimates for theta( 2 ): 0.1171 & ( 0.0149 , 0.3246 )"
## [1] "Point & CI estimates for theta( 3 ): 0.0938 & ( 0.0339 , 0.1815 )"
## [1] "Point & CI estimates for theta( 4 ): 0.1179 & ( 0.0668 , 0.1829 )"
## [1] "Point & CI estimates for theta( 5 ): 0.6068 & ( 0.1612 , 1.3279 )"
## [1] "Point & CI estimates for theta( 6 ): 0.6115 & ( 0.3669 , 0.9046 )"
## [1] "Point & CI estimates for theta( 7 ): 0.8685 & ( 0.0993 , 2.5531 )"
## [1] "Point & CI estimates for theta( 8 ): 0.8635 & ( 0.0982 , 2.4713 )"
## [1] "Point & CI estimates for theta( 9 ): 1.475 & ( 0.4532 , 3.1074 )"
## [1] "Point & CI estimates for theta( 10 ): 1.9471 & ( 1.2253 , 2.8338 )"
labeltext = matrix(c("theta 1", "lower", "mean", "upper","theta 2", "lower", "mean", "upper","theta 3", "lower", "mean", "upper","theta 4", "lower", "mean", "upper","theta 5", "lower", "mean", "upper","theta 6", "lower", "mean", "upper","theta 7", "lower", "mean", "upper","theta 8", "lower", "mean", "upper","theta 9", "lower", "mean", "upper","theta 10", "lower", "mean", "upper"), byrow = TRUE,nrow=10)


forestplot(labeltext,mean = M, lower, upper)

And we can see visually the means and confidence intervals for the failure rates \(\theta_i\). The larger the value of the parameter, the less reliable the pump is of course.

Now assuming we have instead of a fixed \(\alpha\) also a prior on this parameter, such as

\[\alpha \sim\textrm{Exponential}(1.0)\]

Then the posterior becomes

\[h(\underline{\mathbf{\theta}}, \alpha, \beta| \underline{\mathbf{x}}, \underline{\mathbf{t}} ) = \prod_{j=1}^{n=10}{\left[\frac{(t_j\theta_j)^x_j}{x_j!}e^{(-\theta_jt_j)}\frac{\beta^\alpha}{\Gamma(\alpha)}\theta_j^{\alpha-1}e^{(-\beta\theta_j)} \right]} e^{-\alpha}\frac{1}{\Gamma(0.1)}\beta^{-0.9}e^{-\beta}\] The partial condictionals for \(\theta_i\) and \(\beta\) do not change, but for \(\alpha\) it is

\[f_\alpha \propto \frac{\beta^{10\alpha}}{\Gamma(\alpha)^{10}}\prod_{j=1}^{10}\theta_j^{\alpha-1}e^{-\alpha}\]

f.alph = function(alph, bet, thet) {
  f = exp(-a*alph + n*(alph*log(beta)-lgamma(alph))+(alph-1)*sum(log(thet)))
  return(f)
}

# Inits
theta = matrix(0, nrow=N+1, ncol=10) ###### provide for last iteration
alpha[1] = exp(-.4)                      
beta[1] = 1
theta[1,] = rep(0.1, 10)       
accept = 0    ###### change acc to accept
sd = .001     ###### change sd

for(i in 1:N) {
  for(j in 1:10){   
    theta[i+1,j] = rgamma(1, shape = x[j] + alpha[i], rate = t[j] + beta[i])
    ###### update theta[i+1,j]
  }
  beta[i+1] = rgamma(1, shape = b + n*alpha[i],  rate = c + sum(theta[i+1,]))
  ##### update beta
  # M-H for alpha - Method 1  
#  alpha.star = -1
#  while(alpha.star < 0) {
#    alpha.star = rnorm(1,alpha[i],sd)   # sd = 10 may need to tune sd 
#  }
  
  alpha.star=abs(rnorm(1,alpha[i],sd))
  
  f.star = f.alph(alpha.star,beta[i+1],theta[i+1,]) ##### update beta and theta
  fi = f.alph(alpha[i],beta[i+1],theta[i+1,])  ##### update beta and theta
  alpha.prob = min(1,exp(f.star - fi))
  
  if (runif(1) < alpha.prob) {
    alpha[i+1] = alpha.star
    accept = accept+1     ###### change to accept = accept + 1 from accept.alpha = acc + 1
  } else alpha[i+1] = alpha[i]
  #alpha[i+1] = alpha[i]  ##### omit

}
# Acceptance Rate for alpha:

accept/N
## [1] 0.1317
alpha.b = alpha[B:N]
beta.b = beta[B:N]
theta.b = theta[B:N,]

# Parameter Estimates 
mean(alpha.b)
## [1] 1.375494
mean(beta.b)
## [1] 1.878671
apply(theta.b, 2, mean)
##  [1] 0.06593919 0.13684977 0.09891689 0.12008580 0.61547203 0.61385715
##  [7] 0.83926221 0.84275006 1.38163122 1.88948576
apply(theta.b, 2, sd)
##  [1] 0.02615942 0.08906120 0.03917250 0.03070583 0.30497596 0.13719876
##  [7] 0.58251771 0.59493438 0.62625922 0.40386098