/* Change point analysis of coal disaster data using routines from Numerical Recipes */
#include "random.h"
#include <iostream>
#include <iomanip>
#include <fstream>
using namespace std;
ofstream ans("coalnr.ans");
const int n = 112; /* Number of years of data available */
const int m=100; /* Number of replications */
const int t=15; /* Number of iterations */
const int startyear=1851; /* First year for which data is available */
long idum=65537; /* Seed for random number generator */
void graph(int cumY[]); /* Routine to graph cumulative frequency */
static int daytab[2][13] = {
{0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31},
{0, 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31}
};
void gregorian(float a,int &y,int &mo,int &d,
int &h,int &mi,int &s,float &f);
int main (void)
{
int cumY[n],i,j,k,kk,s,sumY;
int year,month,day,hour,minute,second;
double L[n],sumL,p[n],pp[n],sump,cumprob[n+1];
float theta,lambda,b1,b2,U,mean,fraction;
/* Data from B P Carlin, A E Gelfand and A F M Smith, */
/* Hierachical Bayesian Analysis of Changepoint Problems, */
/* Appl. Statist. 41 (1992), 389-405. */
int Y[n]={
4,5,4,1,0,4,3,4,0,6,3,3,4,0,2,6,3,3,5,4,5,3,1,4,4,
1,5,5,3,4,2,5,2,2,3,4,2,1,3,2,1,1,1,1,1,3,0,0,1,0,
1,1,0,0,3,1,0,3,2,2,0,1,1,1,0,1,0,1,0,0,0,2,1,0,0,
0,1,1,0,2,2,3,1,1,2,1,1,1,1,2,4,2,0,0,0,1,4,0,0,0,
1,0,0,0,0,0,1,0,0,1,0,0
};
int a1=1;
int a2=1;
float d1=1;
float d2=1;
sumY=0;
for (i=0; i < n; i++)
{
sumY+=Y[i];
cumY[i]=sumY;
}
graph(cumY);
ans << endl;
for (k=0; k < n; k++)
pp[k]=0;
for (j=0; j < m; j++) /* Replicate m times */
{
k=n/2; /* Initialize k=n/2 */
b1=b2=1; /* Initialize b1=b2=1 */
for (s=1; s < t; s++) /* Iterate t times */
{
/* Sample theta | Y,lambda,b1,b2,k */
theta=gamdev(a1+cumY[k],&idum)/(k+(1/b1));
/* Sample lambda | Y,theta,b1,b2,k */
lambda=gamdev(a2+sumY-cumY[k],&idum)/(n-k+(1/b2));
/* Sample b1 | Y,theta,lambda,b2,k */
b1=(theta+(1/d1))/gamdev(a1,&idum);
/* Sample b2 | Y,theta,lambda,b1,k */
b2=(lambda+(1/d2))/gamdev(a2,&idum);
/* Find L(Y;k,theta,lambda) for k = 0 to n-1 */
sumL=0;
for (k=0; k < n; k++)
{
/* The term -50 is simply to prevent overflow */
L[k]=exp((lambda-theta)*(k+1)+
(log(theta)-log(lambda))*cumY[k]-50);
sumL+=L[k];
}
/* Find p(k | Y,theta,lambda,b1,b2) and cumulation thereof */
sump=0;
cumprob[0]=0;
for (k=0; k < n; k++)
{
p[k]=L[k]/sumL;
sump+=p[k];
cumprob[k+1]=sump;
}
/* Pick U at random between 0 and 1 */
U=ran1(&idum);
/* Sample k | Y,theta,lambda,b1,b2 */
for (i=0; i < n; i++)
if ((cumprob[i] < U)&&(U <= cumprob[i+1])) k=i;
} /* End iteration */
/* Find posterior density and mean of k */
for (i=0; i < n; i++)
pp[i]+=p[i]/m;
} /* End replication */
mean=0;
for (i=0; i < n; i++)
mean+=(startyear+0.5+i)*pp[i];
/* Print out results */
for (i=30; i < 50; i++)
{
ans << startyear+i << " " << pp[i];
ans << endl;
}
ans << endl;
for (i=30; i < 50; i++)
{
ans<< startyear+i << " ";
for (j=1; j < 80; j++)
if (100*pp[i] > j) ans << "*";
ans << endl;
}
ans << endl;
gregorian(mean,year,month,day,hour,minute,second,fraction);
ans << "Mean " << day << "/" << month << "/" << year
<< ", i.e. " << int(mean) << "+" << mean-int(mean)
<< endl;
}
void graph(int cumY[])
{
const int xscale=2;
const int yscale=10;
const int xtick=5;
int xtop,ytop,x,y;
xtop=(n-1)/xscale;
ytop=cumY[n-1]/yscale;
for (y=ytop; y > 0; y--)
{
if (y < 100) ans << " ";
if (y < 10) ans << " ";
ans << 10*y <<" |";
for (x=0; x < xtop; x++)
{
if(y==cumY[xscale*x]/yscale)
ans << "*";
else
ans << " ";
}
ans << endl;
}
ans << " |";
for (x=0; x < xtop; x++)
ans << "_";
ans << endl;
ans << " ";
for (x=0; x < xtop; x++)
if (x==xtick*int((x+1)/xtick)-1)
ans << "|";
else
ans << " ";
ans << endl;
ans << " ";
for (x=0; x < xtop-1; x++)
if (x==5*int(x/5))
ans << startyear+xscale*x+4 << " ";
ans << endl;
}
/* Functions day_of_year and month_day from */
/* B W Kernighan and D M Ritchie, */
/* The C Programming Language, */
/* Englewood Cloffs, NJ: Prentice-Hall 1978, */
/* 1988, Section 5.7. */
/* day_of_year: set day of year from month & day */
int day_of_year(int year, int month, int day)
{
int i, leap;
leap = year%4 == 0 && year%100 != 0 || year%400 == 0;
for (i=1; i < month; i++)
day+=daytab[leap][i];
return day;
}
/* month_day: set month, day from day of year */
void month_day(int year, int yearday, int*pmonth, int *pday)
{
int i, leap;
leap = year%4 == 0 && year%100 != 0 || year%400 == 0;
for (i=1; yearday > daytab[leap][i]; i++)
yearday-=daytab[leap][i];
*pmonth=i;
*pday=yearday;
}
void gregorian(float a,int &y,int &mo,int &d,
int &h,int &mi,int &s,float &f)
{
int yearday,pmonth,pday,leap,daysinyear;
y=int(a);
leap = y%4 == 0 && y%100 != 0 || y%400 == 0;
if (leap)
daysinyear=366;
else
daysinyear=365;
a=daysinyear*(a-int(a));
yearday=int(a);
month_day(y,yearday,&pmonth,&pday);
mo=pmonth;
d=pday;
a=24*(a-yearday);
h=int(a);
a=60*(a-h);
mi=int(a);
a=60*(a-mi);
s=int(a);
f=a-s;
}
/* */