/***************************************************************************
 * Copyright (C) 1995, 1996, Jun Sun, jsun@junsun.net.                     *
 *                                                                         *
 ***************************************************************************
 */
/* implementation of class Binomial */

#include <math.h>
#include "macro.h"
#include "random_variates.h"

Binomial::Binomial(int t, double p)
{ _t = t;
  _p = p;
  bp = new Bernoulli(p);
}

double g(int t,int *f,int n,int M,double mean)
{ int i;
  double sum=0.0;
  double p = mean / (double)t;

  for(i=1;i<=M;i++)
    sum += f[i]*log(t-i+1);
  sum += n*t*log(1-p) + n*mean*log(p/(1-p));
  return sum;
}

Binomial::Binomial(int *data, int num_data)
{ double mean;
  double v;
  double sum,sum2;
  double        max;
  int   *f;
  int   M, M1;
  int  i;

  M = -1;
  sum = 0.0;
  sum2 = 0.0;
  for(i=0;i<num_data;i++)
  { sum += data[i];
    M = MAX(data[i],M);
    sum2 += SQR(data[i]);
  }

  mean = sum / (double) num_data;
  v = sum2/(double)num_data - SQR(mean);

  /* check if these data forms a valid binomial distribution or not */
  if ( mean <= v )
    mythrow("sample data don't form a valid binomial distribution in Binomial::Binomial(int *data, int num_data)");

  /* estimate p and t */
  /* init f */
  f = new int[M+1];
  for(i=0;i<=M;i++) f[i] = 0;
  for(i=0;i<num_data;i++)
    f[data[i]]++;
  for(i=M-1;i>=0;i--)
    f[i] += f[i+1];

  /* compute M1 */
  M1 = FLOOR( mean * (M-1) / (1 - v/mean) );

  /* find _t */
  max = g(M,f,num_data,M,mean);
  _t = M;
  for(i=M+1;i<=M1;i++)
    if (g(i,f,num_data,M,mean) > max)
    { max = g(i,f,num_data,M,mean);
      _t = i;
    }
    else
      break;
  _p = mean / (double)_t;
  bp = new Bernoulli(_p);

  delete[] f;
}


Binomial::~Binomial() 
{ delete bp; }


double Binomial::mass(int x)
{ if ( (x<0) || (x>_t) )
    return 0.0;
  else
    return combination(_t,x)*pow(_p,(double)x)*pow(1-_p,(double)(_t-x));
}


double Binomial::distribution(double x)
{ double f;

  if (x<0.0) return 0;
  if (x>(double)_t) return 1;

  f= 0.0;
  for(int i=0;i<=FLOOR(x);i++)
    f += combination(_t,i)*pow(_p,(double)i)*pow(1-_p,(double)(_t-i));
  return f;
}


int Binomial::sample()
{ int i,sam;
  sam = 0;
  for (i=0;i<_t;i++) sam += bp->sample();
  return sam;
}

int Binomial::binomial_test(int *data, int num_data)
{
  double sum,sum2;
  int   M;
  double mean;
  double variance;
  double v;

  M = -1;
  sum = 0.0;
  sum2 = 0.0;
  for(int i=0;i<num_data;i++)
  { sum += data[i];
    M = MAX(data[i],M);
    sum2 += SQR(data[i]);
  }

  mean = sum / (double) num_data;
  variance = sum2/(double)num_data - SQR(mean);

  /* check if these data forms a valid binomial distribution or not */
  v = (double)(num_data-1)*variance/(double)num_data;
  if ( mean <= v )
    return FALSE;
  else
    return TRUE;
}