# Pumps: Gamma-Poisson Hierarchical Model¶

An example from OpenBUGS [44] and George et al. [36] concerning the number of failures of 10 power plant pumps.

## Model¶

Pump failure are modelled as

where is the number of times that pump failed, and is the operation time of the pump (in 1000s of hours).

## Analysis Program¶

using Distributed
@everywhere using Mamba

## Data
pumps = Dict{Symbol, Any}(
:y => [5, 1, 5, 14, 3, 19, 1, 1, 4, 22],
:t => [94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5]
)
pumps[:N] = length(pumps[:y])

## Model Specification
model = Model(

y = Stochastic(1,
(theta, t, N) ->
UnivariateDistribution[(
lambda = theta[i] * t[i];
Poisson(lambda)) for i in 1:N
],
false
),

theta = Stochastic(1,
(alpha, beta) -> Gamma(alpha, 1 / beta),
true
),

alpha = Stochastic(
() -> Exponential(1.0)
),

beta = Stochastic(
() -> Gamma(0.1, 1.0)
)

)

## Initial Values
inits = [
Dict(:y => pumps[:y], :alpha => 1.0, :beta => 1.0,
:theta => rand(Gamma(1.0, 1.0), pumps[:N])),
Dict(:y => pumps[:y], :alpha => 10.0, :beta => 10.0,
:theta => rand(Gamma(10.0, 10.0), pumps[:N]))
]

## Sampling Scheme
scheme = [Slice([:alpha, :beta], 1.0, Univariate),
Slice(:theta, 1.0, Univariate)]
setsamplers!(model, scheme)

## MCMC Simulations
sim = mcmc(model, pumps, inits, 10000, burnin=2500, thin=2, chains=2)
describe(sim)

## Posterior Predictive Distribution
ppd = predict(sim, :y)
describe(ppd)


## Results¶

## MCMC Simulations

Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2
Samples per chain = 3750

Empirical Posterior Estimates:
Mean         SD        Naive SE        MCSE         ESS
beta 0.93036099 0.517274134 0.00597296721 0.01824153419  804.11618
alpha 0.69679849 0.264420102 0.00305326034 0.00722593007 1339.06428
theta[1] 0.05991674 0.025050913 0.00028926303 0.00032725274 3750.00000
theta[2] 0.10125873 0.078887354 0.00091091270 0.00129985769 3683.18182
theta[3] 0.08910382 0.037307731 0.00043079257 0.00045660986 3750.00000
theta[4] 0.11529009 0.030274265 0.00034957710 0.00032583333 3750.00000
theta[5] 0.59971611 0.316811730 0.00365822675 0.00585119652 2931.65686
theta[6] 0.60967188 0.134761371 0.00155609028 0.00174949908 3750.00000
theta[7] 0.86767451 0.669943846 0.00773584520 0.02858200254  549.40360
theta[8] 0.85445727 0.668132036 0.00771492422 0.02485547216  722.57105
theta[9] 1.55721556 0.753564449 0.00870141275 0.03109274798  587.38464
theta[10] 1.98475207 0.405438843 0.00468160451 0.00912748779 1973.09587

Quantiles:
2.5%       25.0%      50.0%      75.0%      97.5%
beta 0.189558591 0.55034346 0.84013720 1.22683690 2.134324768
alpha 0.286434561 0.50214938 0.66295652 0.85427430 1.319308645
theta[1] 0.021291808 0.04164372 0.05696335 0.07435302 0.118041776
theta[2] 0.008789962 0.04279772 0.08116810 0.13771833 0.305675244
theta[3] 0.032503106 0.06202730 0.08293498 0.11024076 0.176883533
theta[4] 0.063587574 0.09389172 0.11222268 0.13443866 0.182555751
theta[5] 0.153474125 0.36947945 0.54009121 0.77234190 1.364481338
theta[6] 0.371096082 0.51395467 0.60089332 0.69531957 0.898866079
theta[7] 0.077416146 0.37391993 0.71230582 1.17659703 2.616100803
theta[8] 0.072479432 0.35973235 0.69319639 1.16347299 2.655857786
theta[9] 0.463821785 1.01169918 1.43745818 1.96871605 3.351798915
theta[10] 1.269842527 1.70167020 1.95748757 2.24000075 2.861197982

## Posterior Predictive Distribution

Iterations = 2502:10000
Thinning interval = 2
Chains = 1,2
Samples per chain = 3750

Empirical Posterior Estimates:
Mean        SD      Naive SE      MCSE       ESS
y[1]  5.5784000 3.3457654 0.038633571 0.042295418 3750.0000
y[2]  1.6009333 1.7532320 0.020244579 0.022399464 3750.0000
y[3]  5.5854667 3.2976980 0.038078537 0.040503717 3750.0000
y[4] 14.5666667 5.4547523 0.062986054 0.061718435 3750.0000
y[5]  3.1118667 2.4165260 0.027903638 0.035965796 3750.0000
y[6] 19.0780000 6.0320229 0.069651801 0.071700830 3750.0000
y[7]  0.9008000 1.1756177 0.013574864 0.031805581 1366.2355
y[8]  0.8922667 1.1550634 0.013337523 0.025610086 2034.1808
y[9]  3.2585333 2.4039269 0.027758156 0.069046712 1212.1504
y[10] 20.7933333 6.2068800 0.071670876 0.113975236 2965.6896

Quantiles:
2.5% 25.0% 50.0% 75.0% 97.5%
y[1]    1     3     5     8    13
y[2]    0     0     1     2     6
y[3]    1     3     5     7    13
y[4]    5    11    14    18    27
y[5]    0     1     3     4     9
y[6]    9    15    19    23    32
y[7]    0     0     1     1     4
y[8]    0     0     1     1     4
y[9]    0     1     3     5     9
y[10]   10    16    20    25    34