/*

==========================================================================

james.c - a program to run a Markov chain related to James-Stein estimators.
            (Designed for use in conjunction with jpar.c)

Copyright (c) 1999 by Jeffrey S. Rosenthal (jeff@math.toronto.edu).

Licensed for general copying, distribution and modification according to
the GNU General Public License (http://www.gnu.org/copyleft/gpl.html).

----------------------------------------------------

Compile with "cc james.c -lm -o james".

----------------------------------------------------

USAGE:  james [runningtime [seed]]

SYNOPSIS:

james will run for <runningtime> seconds (or 3 seconds if
no runningtime is given), using the drand48() pseudo-random number
generator with intial seed <seed> (or the current time in microseconds
if no seed is given).

During this time, it will continually run a certain Markov chain related
to James-Stein estimators, described in detail in J.S. Rosenthal,
"Analysis of the Gibbs Sampler for a model related to James-Stein
Estimators" (Statistics and Computing, 1996).  It will average the
values of the "theta[0]" variable, but only after discarding the first
139 iterations as burn-in period, and discarding any simulation in
progress to avoid bias.

When done, it will output two numbers: the average of the theta[0]
values considered, and the number of such values averaged.

==========================================================================

*/

#define DEFNUMSECS 3   /* Default number of seconds to run; may be modified. */


#define CONVERGENCETIME 140

#define pi 3.1415926

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <math.h>

main(int argc, char *argv[])
{
    int i, j, numit, runningsecs, tmpseed;
    double tmp;
    double tot = 0.0;
    double N(), IG();
    double mu, A, theta[18], Y[18];
    double Ybar, thetabar;
    int K = 18;
    double V = 0.00434;
    double a = -1;
    double b = 2;
    double thetadev;
    long presenttime();
    long starttime, runningtime;

    /* Determine running time. */
    if ( (argc >= 2) && ( (runningsecs=atoi(argv[1])) > 0 ) )
        ;
    else
        runningsecs = DEFNUMSECS;
    runningtime = runningsecs * 1000000;

    /* Seed the pseudo-random number generator. */
    if ( (argc >= 3) && ( (tmpseed=atoi(argv[2])) > 0 ) )
        srand48(tmpseed);
    else
        seedrand();

    /* Initialise data values Y[i]. */
    Y[0] = 0.395;
    Y[1] = 0.375;
    Y[2] = 0.355;
    Y[3] = 0.334;
    Y[4] = 0.313;
    Y[5] = 0.313;
    Y[6] = 0.291;
    Y[7] = 0.269;
    Y[8] = 0.247;
    Y[9] = 0.247;
    Y[10] = 0.224;
    Y[11] = 0.224;
    Y[12] = 0.224;
    Y[13] = 0.224;
    Y[14] = 0.224;
    Y[15] = 0.200;
    Y[16] = 0.175;
    Y[17] = 0.148;

    /* Compute Ybar. */
    tmp = 0.0;
    for (j=0; j<K; j++)
	tmp = tmp + Y[j];
    Ybar = tmp / K;

    /* Inintialise the theta values. */
    for (j=0; j<K; j++)
	theta[j] = Ybar;

    /* Prepare to begin. */
    starttime = presenttime();
    numit = i = 0;

    /* Loop, adding up Unif[0,1] random variables. */
    while (1) {

	i = i + 1;

	/* Compute thetabar. */
	tmp = 0.0;
	for (j=0; j<K; j++)
	    tmp = tmp + theta[j];
	thetabar = tmp / K;

	/* Compute theta squared-deviation. */
	thetadev = 0.0;
	for (j=0; j<K; j++)
	    thetadev = thetadev + (theta[j]-thetabar)*(theta[j]-thetabar);
	
	/* Run the Markov chain. */
	A = IG( a + 0.5*(K-1), b + 0.5*thetadev );
	mu = N( thetabar, A/K );
	for (j=0; j<K; j++) {
	    theta[j] = N( (mu*V + Y[j]*A) / (V+A), A*V/(V+A) );
	}

	/* Compute average of theta[0]. */
	if (i >= CONVERGENCETIME) {
	    /* Discard simulation in progress if time up. */
	    if ( (numit>0) && (presenttime() > starttime + runningtime) )
		break;
	    tot = tot + theta[0];
	    numit = numit + 1;
	}

    }

    /* Print the result. */
    if (numit > 0)
	printf("%f       %d\n", tot/numit, numit);
    else
	/* Insufficient runs to produce an estimate. */
	printf("0.0      0\n");

}


seedrand()
{
    struct timeval tmptv;
    int seed;

    gettimeofday (&tmptv, (struct timezone *)NULL);
    seed = (int) (tmptv.tv_usec - 1000000 * (int) 
	( ( (double) tmptv.tv_usec ) / ( (double) 1000000.0 ) ) );
    srand48(seed);
    return(0);
}


double IG(double a, double b)
{
    double G();
    return( 1.0 / G(a,b) );
}


double G(double a, double b)
{
    double stdgamma();
    if ( (a<0) || (b<0) ) {
	fprintf(stderr, "G() parameter value out of range.\n");
	exit(1);
    }
    return( stdgamma(a) / b );
}


double stdgamma(double a)  /* gamma with b=1 */
{
    double drand48();
    double b, c, U, V, W, Y, X, Z;
    if (a<1) {
	/* Use Stewart's Theorem, cf. L. Devroye, p. 420, note 3. */
	return ( stdgamma(a+1)*pow(drand48(), 1/a) );
    } else if (a < 1.0001) {
	/* use exponential (a=1) distribution */
	return ( log(1.0/drand48()) );
    } else {
	/* Use Best's rejection algorithm XG, cf. L. Devroye, p. 410. */
        b = a - 1.0;
        c = 3*a - 0.75;
        while(0==0) {
	    U = drand48();
	    V = drand48();
	    W = U*(1-U);
	    Y = sqrt(c/W)*(U-0.5);
	    X = b+Y;
	    if (X >= 0) {
		Z = 64*W*W*W*V*V;
		if ( (Z <= 1 - 2*Y*Y/X) 
			|| (log(Z) <= 2*(b*log(X/b) - Y)) ) {
		    return(X);
		}
	    }
	}
    }
}


double N(double m, double v)
{
    double stdnormvar();
    return( m + sqrt(v)*stdnormvar() );
}


double stdnormvar()
{
    return( sqrt( 2.0*log(1.0/drand48()) ) * cos(2.0*pi*drand48()) );
}


long presenttime()  /* Return present time in microseconds.  */
{
    struct timeval tmp;
    gettimeofday (&tmp, (struct timezone *)NULL);
    return (tmp.tv_usec + 1000000 * tmp.tv_sec);
}