/* Example of rejection sampling (Section 9.5)
 */

#include "random.h"
#include <fstream.h>

const int n=1000;

const int TRUE=1;

const int FALSE=0;

long idum=65537;

ofstream ans("reject.ans");

int main (void)
{
  int i,cont;
  double alpha,beta,c,u,y,x,mean,ss,var,theormean,theorvar;
  alpha=2;
  beta=4;
  c=exp((alpha-1)*log(alpha-1)+(beta-1)*log(beta-1)-
        (alpha+beta-2)*log(alpha+beta-2));
  theormean=alpha/(alpha+beta);
  theorvar=alpha*beta/
           ((alpha+beta)*(alpha+beta)*(alpha+beta+2));
  mean=0;
  ss=0;
  for (i=0; i < n; i++)
  {
     cont=TRUE;
     while (cont)
     {
        y=ran1(&idum);
        u=ran1(&idum);
        if (u <= exp((alpha-1)*log(y)+(beta-1)*log(1-y)))
        {
           x=y;
           mean+=x/n;
           ss+=x*x;
           cont=FALSE;
        }
     }
  }
  var=(ss-n*mean*mean)/(n-1);
  ans << endl;
  ans << " Alpha = " << alpha << " Beta = "     << beta
      << "; Mean = " << mean  << " Variance = " << var << endl;
  ans << " Theoretical values         " << theormean
      << " and        " << theorvar     << endl;
  ans << " Ratios                     " << mean/theormean
      << " and        " << var/theorvar << endl;
  ans << endl;
}

/*