function [lambda, beta] = pumpgibbs(failures, time, alpha, delta, ...
                                    gamma, niter, betastart);
% *************************************************************************
%
% FILE: pumpgibbs.m
%
% Date Created: April 12, 2004
% Author: Mark E. Irwin
%
% Contents: Pump failure example of Gaver and O'Muircheartaigh (1987) as
%           proposed by Gelfand and Smith (1990).  Analysis by Gibbs
%           Sampler.
%
%           p_i | lambda_i ~ Poisson(lambda_i t_i)
%           lambda_i | beta ~ Gamma(alpha, beta)
%           beta | delta, gamma ~ IGamma(gamma, delta)
%
% Revision History
% Date          Name            Changes/Reasons
%
% *************************************************************************
%
% Input:
%
%   failures: number of pump failures
%   time: observation times for the pumps
%   alpha: shape parameter for prior on lambda
%   delta: scale parameter for prior on beta
%   gamma: shape parameter for prior on beta
%   niter: number of Gibbs samples
%   betastart: starting value for beta, the scale parameter for prior on
%              lambda
%
% Output:
%
%   lambda: simulated pump failure rates [niter x npumps]
%   beta: simulated scale parameter for prior on lambda [niter]
%
% *************************************************************************

% Initialize sampler variables

npump = size(failures,1);
ntime = size(time,1);

if (npump ~= ntime) 
    error('Number of pumps does not equal number of times.');
end

lambda = zeros(niter, npump);
beta = zeros(niter,1);

betaold = betastart;

% Run sampler

for iter = 1:niter
    lambdaold = gamrnd(alpha + failures, 1 ./ (time + 1 / betaold));
    lambda(iter, :) = lambdaold';
    betainv = gamrnd(gamma + npump * alpha, 1 / (sum(lambdaold) + delta));
    betaold = 1 / betainv;
    beta(iter) = betaold;
end


