% Title: Football scores an point spreads

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

% ------------------------------------------------------------

% Summary: This is the example of Section 3.4 of Gerlin et al.
% (1995). We have a sample of 672 i.i.d observations from a
% normal distribution with unknown mean and variance. We
% assume a normal prior distribution for the mean and a non-
% informative prior distribution for the variance. The goal
% is to draw a random sample from the joint posterior.

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

% ------------------------------------------------------------

% Variables

% sigma      : variance
% p          : marginal posterior density of the variance
% cdf        : marginal posterior cdf of the variancd
% y          : sample mean
% s          : sample standard error
% mu_0, tau_0: parameters of prior distribution of the mean
% mu_n, tau_n: parameters of posterior conditional distribu-
%              tion of the mean;
% K_1, K_2   : constants
% MEAN       : simulated values for the mean
% VAR        : simulated values for the variance
% i,j        : counters

% Initialization

n = 672;
y = 0.07;
s = 13.86;
tau_0 = 4;
mu_0 = 0;
sigma = 150:.01:250;
cdf = zeros(1,10001);
p = zeros(1,10001);

% Tabulate the distribution of the posterior marginal
% distribution of the variance (using a log transformation).


for i=1:10001;
mu_n=((mu_0/tau_0)+(n*y/sigma(i)))/((1/tau_0)+(n/sigma(i)));
tau_n=1/((1/tau_0)+(n/sigma(i)));
p(i)=0.5*log(tau_n)-(n+2)*0.5*log(sigma(i))-0.5*(((mu_n-mu_0)^2/tau_0)+(((n-1)*s^2+n*(y-mu_n)^2)/sigma(i)));
end;
%Rescaling
K_1=min(p);
p=exp(p-K_1);
K_2=trapz(p);
p=p/K_2;
%Now compute the cdf
cdf=cumsum(p);
% start the iteration
MEAN=zeros(1,1000);
VAR=zeros(1,1000);
for k = 1:1000;
  % draw from the marginal posterior of the variance
  VAR(k) = sigma(min(find(cdf>rand)));
  % draw from the conditional posterior of the mean
  mu_n = ((mu_0/tau_0)+(n*y/VAR(k)))/((1/tau_0) +(n/VAR(k)));
  tau_n = 1/((1/tau_0)+(n/VAR(k)));
  MEAN(k) = mu_n + sqrt(tau_n)*randn;
end;
% End
