htswanalysis/MACS/lib/gsl/gsl-1.11/gsl-randist.c

   1 /* randist/gsl-randist.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 #include <config.h>
  21 #include <stdio.h>
  22 #include <stdlib.h>
  23 #include <math.h>
  24 #include <string.h>
  25
  26 #include <gsl/gsl_randist.h>
  27 #include <gsl/gsl_rng.h>
  28 #include <gsl/gsl_test.h>
  29
  30 void error (const char * s);
  31
  32
  33 int
  34 main (int argc, char *argv[])
  35 {
  36   size_t i,j;
  37   size_t n = 0;
  38   double mu = 0, nu = 0, nu1 = 0, nu2 = 0, sigma = 0, a = 0, b = 0, c = 0;
  39   double zeta = 0, sigmax = 0, sigmay = 0, rho = 0;
  40   double p = 0;
  41   double x = 0, y =0, z=0  ;
  42   unsigned int N = 0, t = 0, n1 = 0, n2 = 0 ;
  43   unsigned long int seed = 0 ;
  44   const char * name ;
  45   gsl_rng * r ;
  46
  47   if (argc < 4)
  48     {
  49       printf (
  50 "Usage: gsl-randist seed n DIST param1 param2 ...\n"
  51 "Generates n samples from the distribution DIST with parameters param1,\n"
  52 "param2, etc. Valid distributions are,\n\n");
  53
  54       printf(
  55 "  beta\n"
  56 "  binomial\n"
  57 "  bivariate-gaussian\n"
  58 "  cauchy\n"
  59 "  chisq\n"
  60 "  dir-2d\n"
  61 "  dir-3d\n"
  62 "  dir-nd\n"
  63 "  erlang\n"
  64 "  exponential\n"
  65 "  exppow\n"
  66 "  fdist\n"
  67 "  flat\n"
  68 "  gamma\n"
  69 "  gaussian-tail\n"
  70 "  gaussian\n"
  71 "  geometric\n"
  72 "  gumbel1\n"
  73 "  gumbel2\n"
  74 "  hypergeometric\n"
  75 "  laplace\n"
  76 "  landau\n"
  77 "  levy\n"
  78 "  levy-skew\n"
  79 "  logarithmic\n"
  80 "  logistic\n"
  81 "  lognormal\n"
  82 "  negative-binomial\n"
  83 "  pareto\n"
  84 "  pascal\n"
  85 "  poisson\n"
  86 "  rayleigh-tail\n"
  87 "  rayleigh\n"
  88 "  tdist\n"
  89 "  ugaussian-tail\n"
  90 "  ugaussian\n"
  91 "  weibull\n") ;
  92       exit (0);
  93     }
  94
  95   argv++ ; seed = atol (argv[0]); argc-- ;
  96   argv++ ; n = atol (argv[0]); argc-- ;
  97   argv++ ; name = argv[0] ; argc-- ; argc-- ;
  98
  99   gsl_rng_env_setup() ;
 100
 101   if (gsl_rng_default_seed != 0) {
 102     fprintf(stderr,
 103             "overriding GSL_RNG_SEED with command line value, seed = %ld\n",
 104             seed) ;
 105   }
 106
 107   gsl_rng_default_seed = seed ;
 108
 109   r = gsl_rng_alloc(gsl_rng_default) ;
 110
 111
 112 #define NAME(x) !strcmp(name,(x))
 113 #define OUTPUT(x) for (i = 0; i < n; i++) { printf("%g\n", (x)) ; }
 114 #define OUTPUT1(a,x) for(i = 0; i < n; i++) { a ; printf("%g\n", x) ; }
 115 #define OUTPUT2(a,x,y) for(i = 0; i < n; i++) { a ; printf("%g %g\n", x, y) ; }
 116 #define OUTPUT3(a,x,y,z) for(i = 0; i < n; i++) { a ; printf("%g %g %g\n", x, y, z) ; }
 117 #define INT_OUTPUT(x) for (i = 0; i < n; i++) { printf("%d\n", (x)) ; }
 118 #define ARGS(x,y) if (argc != x) error(y) ;
 119 #define DBL_ARG(x) if (argc) { x=atof((++argv)[0]);argc--;} else {error( #x);};
 120 #define INT_ARG(x) if (argc) { x=atoi((++argv)[0]);argc--;} else {error( #x);};
 121
 122   if (NAME("bernoulli"))
 123     {
 124       ARGS(1, "p = probability of success");
 125       DBL_ARG(p)
 126       INT_OUTPUT(gsl_ran_bernoulli (r, p));
 127     }
 128   else if (NAME("beta"))
 129     {
 130       ARGS(2, "a,b = shape parameters");
 131       DBL_ARG(a)
 132       DBL_ARG(b)
 133       OUTPUT(gsl_ran_beta (r, a, b));
 134     }
 135   else if (NAME("binomial"))
 136     {
 137       ARGS(2, "p = probability, N = number of trials");
 138       DBL_ARG(p)
 139       INT_ARG(N)
 140       INT_OUTPUT(gsl_ran_binomial (r, p, N));
 141     }
 142   else if (NAME("cauchy"))
 143     {
 144       ARGS(1, "a = scale parameter");
 145       DBL_ARG(a)
 146       OUTPUT(gsl_ran_cauchy (r, a));
 147     }
 148   else if (NAME("chisq"))
 149     {
 150       ARGS(1, "nu = degrees of freedom");
 151       DBL_ARG(nu)
 152       OUTPUT(gsl_ran_chisq (r, nu));
 153     }
 154   else if (NAME("erlang"))
 155     {
 156       ARGS(2, "a = scale parameter, b = order");
 157       DBL_ARG(a)
 158       DBL_ARG(b)
 159       OUTPUT(gsl_ran_erlang (r, a, b));
 160     }
 161   else if (NAME("exponential"))
 162     {
 163       ARGS(1, "mu = mean value");
 164       DBL_ARG(mu) ;
 165       OUTPUT(gsl_ran_exponential (r, mu));
 166     }
 167   else if (NAME("exppow"))
 168     {
 169       ARGS(2, "a = scale parameter, b = power (1=exponential, 2=gaussian)");
 170       DBL_ARG(a) ;
 171       DBL_ARG(b) ;
 172       OUTPUT(gsl_ran_exppow (r, a, b));
 173     }
 174   else if (NAME("fdist"))
 175     {
 176       ARGS(2, "nu1, nu2 = degrees of freedom parameters");
 177       DBL_ARG(nu1) ;
 178       DBL_ARG(nu2) ;
 179       OUTPUT(gsl_ran_fdist (r, nu1, nu2));
 180     }
 181   else if (NAME("flat"))
 182     {
 183       ARGS(2, "a = lower limit, b = upper limit");
 184       DBL_ARG(a) ;
 185       DBL_ARG(b) ;
 186       OUTPUT(gsl_ran_flat (r, a, b));
 187     }
 188   else if (NAME("gamma"))
 189     {
 190       ARGS(2, "a = order, b = scale");
 191       DBL_ARG(a) ;
 192       DBL_ARG(b) ;
 193       OUTPUT(gsl_ran_gamma (r, a, b));
 194     }
 195   else if (NAME("gaussian"))
 196     {
 197       ARGS(1, "sigma = standard deviation");
 198       DBL_ARG(sigma) ;
 199       OUTPUT(gsl_ran_gaussian (r, sigma));
 200     }
 201   else if (NAME("gaussian-tail"))
 202     {
 203       ARGS(2, "a = lower limit, sigma = standard deviation");
 204       DBL_ARG(a) ;
 205       DBL_ARG(sigma) ;
 206       OUTPUT(gsl_ran_gaussian_tail (r, a, sigma));
 207     }
 208   else if (NAME("ugaussian"))
 209     {
 210       ARGS(0, "unit gaussian, no parameters required");
 211       OUTPUT(gsl_ran_ugaussian (r));
 212     }
 213   else if (NAME("ugaussian-tail"))
 214     {
 215       ARGS(1, "a = lower limit");
 216       DBL_ARG(a) ;
 217       OUTPUT(gsl_ran_ugaussian_tail (r, a));
 218     }
 219   else if (NAME("bivariate-gaussian"))
 220     {
 221       ARGS(3, "sigmax = x std.dev., sigmay = y std.dev., rho = correlation");
 222       DBL_ARG(sigmax) ;
 223       DBL_ARG(sigmay) ;
 224       DBL_ARG(rho) ;
 225       OUTPUT2(gsl_ran_bivariate_gaussian (r, sigmax, sigmay, rho, &x, &y),
 226               x, y);
 227     }
 228   else if (NAME("dir-2d"))
 229     {
 230       OUTPUT2(gsl_ran_dir_2d (r, &x, &y), x, y);
 231     }
 232   else if (NAME("dir-3d"))
 233     {
 234       OUTPUT3(gsl_ran_dir_3d (r, &x, &y, &z), x, y, z);
 235     }
 236   else if (NAME("dir-nd"))
 237     {
 238       double *xarr;
 239       ARGS(1, "n1 = number of dimensions of hypersphere");
 240       INT_ARG(n1) ;
 241       xarr = (double *)malloc(n1*sizeof(double));
 242
 243       for(i = 0; i < n; i++) {
 244         gsl_ran_dir_nd (r, n1, xarr) ;
 245         for (j = 0; j < n1; j++) {
 246           if (j) putchar(' ');
 247           printf("%g", xarr[j]) ;
 248         }
 249         putchar('\n');
 250       } ;
 251
 252       free(xarr);
 253     }
 254   else if (NAME("geometric"))
 255     {
 256       ARGS(1, "p = bernoulli trial probability of success");
 257       DBL_ARG(p) ;
 258       INT_OUTPUT(gsl_ran_geometric (r, p));
 259     }
 260   else if (NAME("gumbel1"))
 261     {
 262       ARGS(2, "a = order, b = scale parameter");
 263       DBL_ARG(a) ;
 264       DBL_ARG(b) ;
 265       OUTPUT(gsl_ran_gumbel1 (r, a, b));
 266     }
 267   else if (NAME("gumbel2"))
 268     {
 269       ARGS(2, "a = order, b = scale parameter");
 270       DBL_ARG(a) ;
 271       DBL_ARG(b) ;
 272       OUTPUT(gsl_ran_gumbel2 (r, a, b));
 273     }
 274   else if (NAME("hypergeometric"))
 275     {
 276       ARGS(3, "n1 = tagged population, n2 = untagged population, t = number of trials");
 277       INT_ARG(n1) ;
 278       INT_ARG(n2) ;
 279       INT_ARG(t) ;
 280       INT_OUTPUT(gsl_ran_hypergeometric (r, n1, n2, t));
 281     }
 282   else if (NAME("laplace"))
 283     {
 284       ARGS(1, "a = scale parameter");
 285       DBL_ARG(a) ;
 286       OUTPUT(gsl_ran_laplace (r, a));
 287     }
 288   else if (NAME("landau"))
 289     {
 290       ARGS(0, "no arguments required");
 291       OUTPUT(gsl_ran_landau (r));
 292     }
 293   else if (NAME("levy"))
 294     {
 295       ARGS(2, "c = scale, a = power (1=cauchy, 2=gaussian)");
 296       DBL_ARG(c) ;
 297       DBL_ARG(a) ;
 298       OUTPUT(gsl_ran_levy (r, c, a));
 299     }
 300   else if (NAME("levy-skew"))
 301     {
 302       ARGS(3, "c = scale, a = power (1=cauchy, 2=gaussian), b = skew");
 303       DBL_ARG(c) ;
 304       DBL_ARG(a) ;
 305       DBL_ARG(b) ;
 306       OUTPUT(gsl_ran_levy_skew (r, c, a, b));
 307     }
 308   else if (NAME("logarithmic"))
 309     {
 310       ARGS(1, "p = probability");
 311       DBL_ARG(p) ;
 312       INT_OUTPUT(gsl_ran_logarithmic (r, p));
 313     }
 314   else if (NAME("logistic"))
 315     {
 316       ARGS(1, "a = scale parameter");
 317       DBL_ARG(a) ;
 318       OUTPUT(gsl_ran_logistic (r, a));
 319     }
 320   else if (NAME("lognormal"))
 321     {
 322       ARGS(2, "zeta = location parameter, sigma = scale parameter");
 323       DBL_ARG(zeta) ;
 324       DBL_ARG(sigma) ;
 325       OUTPUT(gsl_ran_lognormal (r, zeta, sigma));
 326     }
 327   else if (NAME("negative-binomial"))
 328     {
 329       ARGS(2, "p = probability, a = order");
 330       DBL_ARG(p) ;
 331       DBL_ARG(a) ;
 332       INT_OUTPUT(gsl_ran_negative_binomial (r, p, a));
 333     }
 334   else if (NAME("pareto"))
 335     {
 336       ARGS(2, "a = power, b = scale parameter");
 337       DBL_ARG(a) ;
 338       DBL_ARG(b) ;
 339       OUTPUT(gsl_ran_pareto (r, a, b));
 340     }
 341   else if (NAME("pascal"))
 342     {
 343       ARGS(2, "p = probability, n = order (integer)");
 344       DBL_ARG(p) ;
 345       INT_ARG(N) ;
 346       INT_OUTPUT(gsl_ran_pascal (r, p, N));
 347     }
 348   else if (NAME("poisson"))
 349     {
 350       ARGS(1, "mu = scale parameter");
 351       DBL_ARG(mu) ;
 352       INT_OUTPUT(gsl_ran_poisson (r, mu));
 353     }
 354   else if (NAME("rayleigh"))
 355     {
 356       ARGS(1, "sigma = scale parameter");
 357       DBL_ARG(sigma) ;
 358       OUTPUT(gsl_ran_rayleigh (r, sigma));
 359     }
 360   else if (NAME("rayleigh-tail"))
 361     {
 362       ARGS(2, "a = lower limit, sigma = scale parameter");
 363       DBL_ARG(a) ;
 364       DBL_ARG(sigma) ;
 365       OUTPUT(gsl_ran_rayleigh_tail (r, a, sigma));
 366     }
 367   else if (NAME("tdist"))
 368     {
 369       ARGS(1, "nu = degrees of freedom");
 370       DBL_ARG(nu) ;
 371       OUTPUT(gsl_ran_tdist (r, nu));
 372     }
 373   else if (NAME("weibull"))
 374     {
 375       ARGS(2, "a = scale parameter, b = exponent");
 376       DBL_ARG(a) ;
 377       DBL_ARG(b) ;
 378       OUTPUT(gsl_ran_weibull (r, a, b));
 379     }
 380   else
 381     {
 382       fprintf(stderr,"Error: unrecognized distribution: %s\n", name) ;
 383     }
 384
 385   return 0 ;
 386 }
 387
 388
 389 void
 390 error (const char * s)
 391 {
 392   fprintf(stderr, "Error: arguments should be %s\n",s) ;
 393   exit (EXIT_FAILURE) ;
 394 }