/* Hierachical normal model at end of Section 9.2
 */

#include "random.h"
#include <iostream>
#include <iomanip>
#include <fstream.h>
using namespace std;

const int niter=25;

const int r = 4;
const int n[4]={4,6,6,8};

double x[4][8] = {
  {62.0, 60.0, 63.0, 59.0,  0.0,  0.0,  0.0,  0.0},
  {63.0, 67.0, 71.0, 64.0, 65.0, 66.0,  0.0,  0.0},
  {68.0, 66.0, 71.0, 67.0, 68.0, 68.0,  0.0,  0.0},
  {56.0, 62.0, 60.0, 61.0, 63.0, 64.0, 63.0, 59.0}
};

char pagethrow=12;

int main(void)
{
   ofstream ans("hiernorm.ans");
   ans.setf(ios::fixed);
   ans << endl;
   ans << "Based on A. Gelman, J.B. Carlin, "
       << "H.S. Stern and D.B. Rubin"<< endl;
   ans << "Bayesian Data Analysis, London: Chapman & Hall "
       << "1995, Sec. 9.2" << endl;
   int i,j,t,N;
   double xidot[r],ssi[r],theta[r],v[r];
   double xdotdot,ssw,ssb;
   double mu,phi,psi,muold,phiold,psiold;
   N=0;
   for (i = 0; i < r; i++)
   {
      N+=n[i];
      xidot[i] = 0;
      ssi[i] = 0;
      for (j = 0; j < n[i]; j++)
         xidot[i]+=x[i][j]/n[i];
   }
   xdotdot=0;
   for (i = 0; i < r; i++)
      for (j = 0; j < n[i]; j++)
      {
         xdotdot+=x[i][j]/N;
         ssi[i]+=(x[i][j]-xidot[i])*(x[i][j]-xidot[i]);
      }
   ssw=ssb=0;
   for (i = 0; i < r; i++)
   {
      ssw+=ssi[i];
      ssb+=(xidot[i]-xdotdot)*(xidot[i]-xdotdot);
   }
   mu=muold=xdotdot;
   phi=phiold=ssw/(N-1);
   psi=psiold=ssb/(r-1);
   for (t = 0; t < niter; t++)
   {
      muold=mu;
      phiold=phi;
      psiold=psi;
      mu=phi=psi=0;
      for (i = 0; i < r; i++)
      {
         v[i]=1/(1/psiold+n[i]/phiold);
         theta[i]=v[i]*(muold/psiold+n[i]*xidot[i]/phiold);
      }
      mu=0;
      for (i = 0; i < r; i++)
         mu+=theta[i]/r;
      for (i = 0; i < r; i++)
         for (j = 0; j < n[i]; j++)
            phi+=(v[i]+(x[i][j]-theta[i])*
                 (x[i][j]-theta[i]))/(N+2);
      for (i = 0; i < r; i++)
         psi+=(v[i]+(mu-theta[i])*(mu-theta[i]))/r;
   }
   ans << endl;
   for (i = 0; i < r; i++)
      ans << "Theta["<