% This is the example of Section 3.2 of Gelman et al. (1995):
% Estimating the Speed of Light

% Author: Patrick Waelbroeck, P.O. Box 202376, New Haven, CT
% 06520, USA; pwaelbro@pantheon.yale.edu

% The sample mean is 26.21, while the sample standard
% deviation is 10.75; there are 66 observations.

% Reference

% Gelman, A., Carlin, J.B., Stern, H.S., and Rubin, D.B.
% (1995), Bayesian Data Analysis, Chapman & Hall.
% ------------------------------------------------------------

% Begin

% External routine: randgam(n1, n2, alpha, beta)

% Variables

% n  : counter (1000 simulated values)
% m  : simulated posterior distribution of the mean
% s  : simulated posterior distribution of the variance
% m1 : to compute the central posterior confidence interval

% Initialization

s = zeros(1000,1);
m = zeros(1000,1);

% Start the loop

for n=1:1000;

  % Random draw from inverse chi-squared(65, 115.56)
  s(n) = 1 / randgam(1,1,32.5,3755.7);

  % Random draw from normal(26.21, s(n)/66)
  m(n) = 26.21 + sqrt(s(n)/66)*randn;

end; % End for n

% Compute the 95% central posterior confidence interval

m1 = sort(m);

fprintf('Confidence interval: [%.2f; %.2f]\n', m1(25), m1(976));

% End
