Re: Error in gsl_ran_gaussian_tail ?!

Re: Error in gsl_ran_gaussian_tail ?!

[Brian Gough]

Achim Gaedke writes:

 > 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: 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.).

Thanks, the case a<0 was not implemented in randist/gausstail.c.

It would make sense for it to always return the upper tail, which we
can easily do by removing the 'fabs' from the first branch of the code
there.  I will do that.

bjg|debian> cvs diff -u gausstail.c
Index: gausstail.c
===================================================================
RCS file: /cvs/gsl/gsl/randist/gausstail.c,v
retrieving revision 1.4
diff -u -r1.4 gausstail.c
--- gausstail.c 2001/04/23 09:38:28     1.4
+++ gausstail.c 2001/06/23 11:33:35
@@ -27,9 +27,8 @@
 double
 gsl_ran_gaussian_tail (const gsl_rng * r, const double a, const double sigma)
 {
-  /* Returns a gaussian random variable larger than s
-   * This implementation does one-tailed deviates.
-   * FIXME: what to do about a < 0?
+  /* Returns a gaussian random variable larger than a
+   * This implementation does one-sided upper-tailed deviates.
    */
 
   double s = a / sigma;
@@ -43,7 +42,7 @@
 
       do
        {
-         x = fabs (gsl_ran_gaussian (r, 1.0));
+         x = gsl_ran_gaussian (r, 1.0);
        }
       while (x < s);
       return x * sigma;

Error in gsl_ran_gaussian_tail ?!

[Achim Gaedke]

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #1: Type: text/plain, Size: 1729 bytes --]

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;
    }
}