/* Generate a number of Poission-distrubuted data, and
   compare the straight average with one where weights 
   are chosen as 1/sigma^2, with sigma=sqrt(N) */
/* compile with gcc -lgsl -lgslcblas testpoisson.c */
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
int main (void)
{
  int ntrial=100,nsample=100;
  double mu=100.;
  int i,j;
  double f,s,w,sumf,sumwf,sumw,av,wav;
  /* Random number generator setup; see example 
http://www.gnu.org/software/gsl/manual/html_node/Random-Number-Distribution-Examples.html#Random-Number-Distribution-Examples
  */
  const gsl_rng_type * T;
  gsl_rng * r;
  gsl_rng_env_setup();
  T = gsl_rng_default;
  r = gsl_rng_alloc (T);
  /* now the real action */
  /* one hundred times, ... */
  for(i=0;i<ntrial;i++){
    /* average 100 samples */
    sumf = sumw = sumwf = 0;
    for(j=0;j<nsample;j++){
      f = (double) gsl_ran_poisson(r, mu);
      s = sqrt(f);
      w = pow((1./s),2);
      sumf = sumf+f;
      sumw = sumw+w;
      sumwf = sumwf+w*f;
    }
    av = sumf/(double)nsample;
    wav = sumwf/sumw;
    printf("%3d %10.6f %10.6f\n",i, av, wav);
  }
}