/* Genetic linkage example from Sections 9.2 and 9.3
 */

#include "random.h"
#include <iostream>
#include <iomanip>
#include <fstream>
using namespace.std

const int niter = 50;

const int minitial = 20;

int y[3000];

const int etadivs = 1000;

const int mag = 10;

const int scale = 8;

const int  x[5] = {0,125,18,20,34};

long idum = 65537;

char pagethrow=12;

int main(void)
{
   ofstream ans("dataaug.ans");
   ans.setf(ios::fixed);
   int i,i0,j,t,k,m,mold,etastep;
   double u,v,eta;
   double p[etadivs+1];
   m = minitial;
   for (i = 0; i < m; i++)
   {
      eta=ran1(&idum);
      y[i]=int(bnldev(eta/(eta+2),x[1],&idum));
   }
   for (t = 1; t <= niter; t++)
   {
      mold = m;
      if (t > 30)
         m = 200;
      if (t > 40)
         m = 1000;
      for (i = 0; i < m; i++)
      {
         i0=int(mold*ran1(&idum));
         u=gamdev(y[i0]+x[4]+1,&idum);
         v=gamdev(x[2]+x[3]+1,&idum);
         eta=u/(u+v);
         y[i]=int(bnldev(eta/(eta+2),x[1],&idum));
      }
   }
   for (etastep = 0; etastep <= etadivs; etastep++)
   {
      eta=float(etastep)/float(etadivs);
      p[etastep]=0;
      for (i = 0; i < m; i++)
      {
         p[etastep]+=exp((y[i]+x[4])*log(eta)+
                         (x[2]+x[3])*log(1-eta)+
                         gammln(y[i]+x[2]+x[3]+x[4]+2)-
                         gammln(y[i]+x[4]+1)-
                         gammln(x[2]+x[3]+1))/m;
      }
      if ((0.62 <= eta)&&(eta <= 0.635))
         ans << " p["  << setprecision(3) << eta
             << "] = " << setprecision(4) << p[etastep] << endl;
   }
   ans << setprecision(1) << endl;
   for (etastep = 0; etastep <= etadivs; etastep+=mag)
   {
      eta=float(etastep)/float(etadivs);
      if ((0.45 <= eta)&&(eta <= 0.81))
      {
          if (10*eta-int(10*eta+0.001) < 0.001) ans << eta;
          else ans << "   ";
          ans << "|";
          for (k = 1; k <= scale*p[etastep]; k++) ans << "*";
          ans << endl;
      }
   }
}

/*