This is the mail archive of the gsl-discuss@sources.redhat.com mailing list for the GSL project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]

Error in gsl_ran_gaussian_tail ?!


Hello!

I've tried to generate random numbers distributed according to a gaussian
tail distribution...

I've set the cutoff a to -4, 0 and 4. For -4 and 0 I recieved the same
result, the pictures are made with gnuplot:

 http://www.zpr.Uni-Koeln.DE/~achim/gsl_ran_test/gsl_test_0_7.png

This seems to be 2*gauss at positve values  and is correct for cutoff 0.
cutoff 4 is correct.
I just need negative cutoffs (random money deposits e.g.).

The source code (appended) and pictures are available at:

http://www.zpr.Uni-Koeln.DE/~achim/gsl_ran_test

I've ran tests on Redhat 7.1 Linux both on gsl-0.7 and gsl-0.8.

Thanks for any hints!

Achim

Achim Gaedke, ZPR
Weyertal 80, 50931 Köln
Tel: +49 221 470 6021

#include <stdio.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_histogram.h>


int main(void)
{
  FILE* out_file; /* output stream*/
  int i;
  double cut_off;
  double sigma;
  const int sample_no=100000;   /* sample numbers */
  const int bin_no=100;      /* number of bins in histogram */
  gsl_histogram* my_hist;
  gsl_rng* my_rng;

  out_file=stdout;
  my_hist=gsl_histogram_calloc_uniform(bin_no,-10,10);
  my_rng=gsl_rng_alloc(gsl_rng_default);

  sigma=1;
  cut_off=-4;

  fprintf(out_file,"plot '-' using 1:3 title '1st',");
  fprintf(out_file,"'-' using 1:3 title '2nd',");
  fprintf(out_file,"'-' using 1:3 title '3rd'\n");
  while (cut_off<=4)
    {
      fprintf(out_file,"# cut_off %f, sigma %f\n",cut_off,sigma);
      gsl_histogram_reset(my_hist);
      for (i=0;i<sample_no;i++)
	gsl_histogram_increment(my_hist,gsl_ran_gaussian_tail(my_rng,cut_off,sigma));
      gsl_histogram_fprintf(out_file,my_hist,"%f","%f");
      fprintf(out_file,"e\n");
      cut_off+=4;
    }
}


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]