/***************************************************************************
 * Copyright (C) 1995, 1996, Jun Sun, jsun@junsun.net.                     *
 *                                                                         *
 ***************************************************************************
 */
/* implementation of negtive poisson distribution */

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


Poisson::Poisson(double lamda)
{ _lamda = lamda;
}

Poisson::Poisson(int *data, int num_data)
{ double sum;

  sum = 0.0;
  for(int i=0;i<num_data;i++)
    sum += data[i];

  _lamda = sum / (double) num_data;
}

Poisson::~Poisson() {}

double Poisson::mass(int x)
{ if (x<0)
    return 0.0;
  else
    return exp(-_lamda) * pow(_lamda,(double)x) / factorial(x);
}


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

  if (x<0.0) return 0.0;

  f= 0.0;
  for(int i=0;i<=FLOOR(x);i++)
    f += pow(_lamda,(double)i) / factorial(i);
  return exp(-_lamda) * f;
}


int Poisson::sample()
{ int i=0;
  double a = exp(-_lamda);
  double b = 1.0;

  do
  {
    b *= randp->sample();
    i++;
  } while(b>=a);
  
  return i-1;
}