/* Genetic linkage example from Chapter 9, Exercise 3
 */

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

const int niter = 50;

const int minitial = 20;

int y[5][3000];

const int etadivs = 1000;

const int mag = 10;

const int scale = 4;

const int x[5] = {0,461,130,161,515};

long idum = 65537;

char pagethrow=12;

int main(void)
{
   ofstream ans("dataaug2.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[1][i]=int(bnldev(eta/(eta+1),x[1],&idum));
      y[4][i]=int(bnldev((1-eta)/(2-eta),x[4],&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[1][i0]+x[2]+1,&idum);
         v=gamdev(x[3]+y[4][i0]+1,&idum);
         eta=u/(u+v);
         y[1][i]=int(bnldev(eta/(eta+1),x[1],&idum));
         y[4][i]=int(bnldev((1-eta)/(2-eta),x[4],&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[1][i]+x[2])*log(eta)+
                         (x[3]+y[4][i])*log(1-eta)+
                         gammln(y[1][i]+x[2]+x[3]+y[4][i]+2)-
                         gammln(y[1][i]+x[2]+1)-
                         gammln(x[3]+y[4][i]+1))/m;
      }
      if ((0.43 <= eta)&&(eta <= 0.45))
         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.3 <= eta)&&(eta <= 0.55))
      {
          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;
      }
   }
}

/*