/* Chained data augmentation - Example from Casella and George
*/
#include "random.h"
#include <iostream>
#include <iomanip>
#include <fstream>
#include <cmath>
using namespace std;
const int nr = 50;
const int m = 500;
const int k = 10;
const int n = 16;
const int alpha = 2;
const int beta = 4;
const int lambda = 16;
const int maxn = 24;
char pagethrow=12;
double u,v,y,term,ratio,compx2,x2h,x2f;
int ifail,i,j,t,x,newalpha,newbeta,diff,nn,practmaxn;
long idum=65537;
int h[n+1],bb[n+1],f[n+1];
int hp[maxn+1],fp[maxn+1];
double bbe[n+1],fe[n+1];
double fep[maxn+1];
double betabinomial(int x, int n, double alpha, double beta)
{
y = log(bico(n,x));
y = y + gammln(alpha + beta) - gammln(alpha) - gammln(beta);
y = y + gammln(x + alpha) + gammln(n - x + beta) -
gammln(alpha + beta + n);
y = exp(y);
return y;
}
int main(void)
{
ofstream gibbsans("chained.ans");
gibbsans.setf(ios::fixed);
gibbsans << "Based on 'Explaining the Gibbs sampler', ";
gibbsans << "C. Casella" << endl;
gibbsans << "and E.I. George, Amer. Statist. ";
gibbsans << "46 (3) (1992), 167-174." << endl;
gibbsans << endl;
ifail = 0;
for (t = 0; t < n+1; t++)
{
h[t] = 0;
fe[t] = 0;
}
for (i = 0; i < m; i++)
{
y = ran1(&idum);
for (j = 0; j < k; j++)
{
x = int(bnldev(y,n,&idum));
newalpha = x + alpha;
newbeta = n - x + beta;
u = gamdev(newalpha,&idum);
v = gamdev(newbeta,&idum);
y = u/(u + v);
}
for (t = 0; t < n+1; t++)
{
if (t == x)
h[t] = h[t] + 1;
term = bico(n,t)*exp(t*log(y)+(n-t)*log(1-y));
fe[t] = fe[t] + term;
}
}
gibbsans << "Histogram (cf. Fig. 1))" << endl;
gibbsans << endl;
gibbsans << " t Obs Exp Diff Ratio ";
gibbsans << "Comp of X2" << endl;
gibbsans << endl;
x2h = 0;
for (t = 0; t < n + 1; t++)
{
bbe[t] = m*betabinomial(t,n,alpha,beta);
bb[t] = int(bbe[t]);
diff = h[t] - bb[t];
ratio = h[t]/bbe[t];
compx2 = (h[t]-bbe[t])*(h[t]-bbe[t])/bbe[t];
x2h = x2h + compx2;
gibbsans << setw(10) << t << setw(10) << h[t]
<< setw(10) << bb[t] << setw(10) << diff;
gibbsans << setw(11) << setprecision(3) << ratio;
gibbsans << setw(13) << setprecision(3) << compx2;
gibbsans << endl;
}
gibbsans << endl;
gibbsans << "Chi-squared equals ";
gibbsans << x2h;
gibbsans << " on " << n << " degrees of freedom" << endl;
gibbsans << endl;
gibbsans << "Estimated densities (cf. Fig. 3)" << endl;
gibbsans << endl;
gibbsans << " t Obs Exp Diff Ratio ";
gibbsans << "Comp of X2" << endl;
gibbsans << endl;
x2f = 0;
for (t = 0; t < n+1; t++)
{
f[t] = int(fe[t]);
diff = f[t] - bb[t];
ratio = f[t]/bbe[t];
compx2 = (f[t]-bbe[t])*(f[t]-bbe[t])/bbe[t];
x2f = x2f + compx2;
gibbsans << setw(10) << t << setw(10)
<< f[t] << setw(10)
<< setw(10) << bb[t] << setw(10) << diff;
gibbsans << setw(11) << setprecision(3) << ratio;
gibbsans << setw(13) << setprecision(3) << compx2
<< endl;
}
gibbsans << endl;
gibbsans << "Chi-squared equals ";
gibbsans << x2f;
gibbsans << " on " << n << " degrees of freedom." << endl;
for (t = 0; t < maxn+1; t++)
{
hp[t] = 0;
fep[t] = 0;
}
for (i = 0; i < m; i++)
{
y = 0.5;
nn = int((1-y)*lambda);
for (j = 0; j < k; j++)
{
x = int(bnldev(y,nn,&idum));
newalpha = x + alpha;
newbeta = nn - x + beta;
u = gamdev(newalpha,&idum);
v = gamdev(newbeta,&idum);
y = u/(u + v);
nn = x + int(poidev((1-y)*lambda,&idum));
}
for (t = 0; t < maxn+1; t++)
{
if (t == x)
hp[t] = hp[t] + 1;
if (t <= nn)
{
term = bico(nn,t)*exp(t*log(y)+
(nn-t)*log(1-y));
fep[t] = fep[t] + term;
}
}
}
gibbsans << pagethrow;
gibbsans << "Histogram (n random)" << endl;
gibbsans << endl;
gibbsans << " t Obs Histogram" << endl;
gibbsans << endl;
practmaxn = 4*n/3;
for (t = 0; t < practmaxn+1; t++)
{
gibbsans << setw(2) << t << setw(10) << hp[t]
<< " ";
for (j = 0; j < hp[t]; j++) gibbsans << "*";
gibbsans << endl;
}
gibbsans << endl;
gibbsans << "Estimated densities (n random; cf. Fig. 5)"
<< endl;
gibbsans << endl;
gibbsans << " t Obs Estimate" << endl;
gibbsans << endl;
x2f = 0;
for (t = 0; t < practmaxn+1; t++)
{
fp[t] = int(fep[t]);
gibbsans << setw(2) << t << setw(10) << fp[t]
<< " ";
for (j = 0; j < fp[t]; j++) gibbsans << "*";
gibbsans << endl;
}
}
/* */