/* Rat data in Chapter 9, Exercise 11 */
#include "random.h"
#include "jacobi.h"
#include "wishart.h"
#include <iostream.h>
#include <iomanip.h>
#include <fstream.h>
using namespace std;
#define P 2
int p=P;
#include "wishart.c" /* Must follow #define P */
const int FALSE=0;
const int TRUE=1;
const int m = 500;
const int k = 30;
const int ni = 5;
long idum=65537;
float epsilon = 0.001;
int main (void)
{
int i,ii,j;
float y[k][ni];
float x[] = {8, 15, 22, 29, 36};
float xbar=0.0;
float sum_x2=0.0;
float a,aalpha,abeta,b,balpha,bbeta,
alpha0,beta0,alphabar,betabar;
float alpha[k],beta[k];
float sigma2,sigmaa2,sigmab2;
float var[P][P],w[P][P],invw[P][P],invsigma[P][P];
float vara,varb,suma,sumb,meana,meanb,z;
for (j=0; j < ni; j++)
{
xbar+=x[j]/ni;
sum_x2+=x[j]*x[j];
}
sum_x2-=ni*xbar*xbar;
ofstream ans("rats.ans");
ifstream data("rats.dat");
for (i=0; i < k; i++)
{
data >> ii;
for (j=0; j < ni; j++)
data >> y[i][j];
}
alpha0=beta0=0.0;
a=aalpha=abeta=epsilon; /* B P Carlin and T A Louis p. 169 */
b=balpha=bbeta=1/epsilon; /* B P Carlin and T A Louis p. 170 */
sigma2=1.0; /* Initially sigma2 is IG(a,b) */
sigmaa2=100; /* A E Gelfand et al. p. 979 col.1 */
sigmab2=0.1; /* A E Gelfand et al. p. 979 col.1 */
/* Thus R = (100 0 ) */
/* ( 0 0.1) */
/* Take values for alpha[i] and beta[i] given */
/* alpha0, beta0, sigmaa2, sigmab2 and sigma2 */
alphabar=0.0;
betabar=0.0;
for (i=0; i < k; i++)
{
vara=ni/sigma2 + 1/sigmaa2;
suma=0.0;
for (j=0; j < ni; j++) suma+=y[i][j];
meana=(suma/sigma2 + alpha0/sigmaa2)/vara;
/* alpha[i] ~ N(meana,vara) */
alpha[i]=meana+sqrt(vara)*gasdev(&idum);
alphabar+=alpha[i]/k;
varb=sum_x2/sigma2 + 1/sigmab2;
sumb=0.0;
for (j=0; j < ni; j++) sumb+=(x[j]-xbar)*y[i][j];
meanb=(sumb/sigma2 + 1/sigmab2)/varb;
/* beta[i] ~ N(meanb,varb) */
beta[i]=meanb+sqrt(varb)*gasdev(&idum);
betabar+=beta[i]/k;
}
/* Initialize var (capital sigma) */
var[0][0]=sigmaa2;
var[0][1]=0;
var[1][0]=0;
var[1][1]=sigmab2;
/* Take values for alpha0 and beta0 given */
/* alpha[i], beta[i], sigmaa2, sigmab2 and sigma2 */
wish(p,k,var,w,invw,&idum);
/* alpha0 ~ N(alphabar,sigmaa2/k) */
alpha0=alphabar+sqrt(sigmaa2/k)*gasdev(&idum);
/* beta0 ~ N(betabar,sigmab2/k) */
beta0=betabar+sqrt(sigmab2/k)*gasdev(&idum);
/* See last displayed formula on p. 168 of */
/* B P Carlin and T A Louis */
/* simplified by taking C^{-1} = 0 */
/* Take values for sigmaa2 and sigmab2 given */
/* alpha[i], beta[i], alpha0, beta0 and sigma2 */
/* sigmaa2 ~ IG(alpha0,beta0) */
sigmaa2=1/(beta0*gammadev(alpha0,&idum));
/* sigmab2 ~ IG(alpha0,beta0) */
sigmab2=1/(beta0*gammadev(alpha0,&idum));
/* Take value for sigma2 given */
/* alpha[i], beta[i], alpha0, beta0, sigmaa2 and sigmab2 */
/* sigma2 ~ IG(alpha0,beta0) */
sigma2=1/(beta0*gammadev(alpha0,&idum));
}
/* */