diff options
Diffstat (limited to 'gsl-1.9/gsl-randist.c')
-rw-r--r-- | gsl-1.9/gsl-randist.c | 393 |
1 files changed, 393 insertions, 0 deletions
diff --git a/gsl-1.9/gsl-randist.c b/gsl-1.9/gsl-randist.c new file mode 100644 index 0000000..714493d --- /dev/null +++ b/gsl-1.9/gsl-randist.c @@ -0,0 +1,393 @@ +/* randist/gsl-randist.c + * + * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <stdio.h> +#include <stdlib.h> +#include <math.h> +#include <string.h> + +#include <gsl/gsl_randist.h> +#include <gsl/gsl_rng.h> +#include <gsl/gsl_test.h> + +void error (const char * s); + + +int +main (int argc, char *argv[]) +{ + size_t i,j; + size_t n = 0; + double mu = 0, nu = 0, nu1 = 0, nu2 = 0, sigma = 0, a = 0, b = 0, c = 0; + double zeta = 0, sigmax = 0, sigmay = 0, rho = 0; + double p = 0; + double x = 0, y =0, z=0 ; + unsigned int N = 0, t = 0, n1 = 0, n2 = 0 ; + unsigned long int seed = 0 ; + const char * name ; + gsl_rng * r ; + + if (argc < 4) + { + printf ( +"Usage: gsl-randist seed n DIST param1 param2 ...\n" +"Generates n samples from the distribution DIST with parameters param1,\n" +"param2, etc. Valid distributions are,\n" +"\n" +" beta\n" +" binomial\n" +" bivariate-gaussian\n" +" cauchy\n" +" chisq\n" +" dir-2d\n" +" dir-3d\n" +" dir-nd\n" +" erlang\n" +" exponential\n" +" exppow\n" +" fdist\n" +" flat\n" +" gamma\n" +" gaussian-tail\n" +" gaussian\n" +" geometric\n" +" gumbel1\n" +" gumbel2\n" +" hypergeometric\n" +" laplace\n" +" landau\n" +" levy\n" +" levy-skew\n" +" logarithmic\n" +" logistic\n" +" lognormal\n" +" negative-binomial\n" +" pareto\n" +" pascal\n" +" poisson\n" +" rayleigh-tail\n" +" rayleigh\n" +" tdist\n" +" ugaussian-tail\n" +" ugaussian\n" +" weibull\n") ; + exit (0); + } + + argv++ ; seed = atol (argv[0]); argc-- ; + argv++ ; n = atol (argv[0]); argc-- ; + argv++ ; name = argv[0] ; argc-- ; argc-- ; + + gsl_rng_env_setup() ; + + if (gsl_rng_default_seed != 0) { + fprintf(stderr, + "overriding GSL_RNG_SEED with command line value, seed = %ld\n", + seed) ; + } + + gsl_rng_default_seed = seed ; + + r = gsl_rng_alloc(gsl_rng_default) ; + + +#define NAME(x) !strcmp(name,(x)) +#define OUTPUT(x) for (i = 0; i < n; i++) { printf("%g\n", (x)) ; } +#define OUTPUT1(a,x) for(i = 0; i < n; i++) { a ; printf("%g\n", x) ; } +#define OUTPUT2(a,x,y) for(i = 0; i < n; i++) { a ; printf("%g %g\n", x, y) ; } +#define OUTPUT3(a,x,y,z) for(i = 0; i < n; i++) { a ; printf("%g %g %g\n", x, y, z) ; } +#define INT_OUTPUT(x) for (i = 0; i < n; i++) { printf("%d\n", (x)) ; } +#define ARGS(x,y) if (argc != x) error(y) ; +#define DBL_ARG(x) if (argc) { x=atof((++argv)[0]);argc--;} else {error( #x);}; +#define INT_ARG(x) if (argc) { x=atoi((++argv)[0]);argc--;} else {error( #x);}; + + if (NAME("bernoulli")) + { + ARGS(1, "p = probability of success"); + DBL_ARG(p) + INT_OUTPUT(gsl_ran_bernoulli (r, p)); + } + else if (NAME("beta")) + { + ARGS(2, "a,b = shape parameters"); + DBL_ARG(a) + DBL_ARG(b) + OUTPUT(gsl_ran_beta (r, a, b)); + } + else if (NAME("binomial")) + { + ARGS(2, "p = probability, N = number of trials"); + DBL_ARG(p) + INT_ARG(N) + INT_OUTPUT(gsl_ran_binomial (r, p, N)); + } + else if (NAME("cauchy")) + { + ARGS(1, "a = scale parameter"); + DBL_ARG(a) + OUTPUT(gsl_ran_cauchy (r, a)); + } + else if (NAME("chisq")) + { + ARGS(1, "nu = degrees of freedom"); + DBL_ARG(nu) + OUTPUT(gsl_ran_chisq (r, nu)); + } + else if (NAME("erlang")) + { + ARGS(2, "a = scale parameter, b = order"); + DBL_ARG(a) + DBL_ARG(b) + OUTPUT(gsl_ran_erlang (r, a, b)); + } + else if (NAME("exponential")) + { + ARGS(1, "mu = mean value"); + DBL_ARG(mu) ; + OUTPUT(gsl_ran_exponential (r, mu)); + } + else if (NAME("exppow")) + { + ARGS(2, "a = scale parameter, b = power (1=exponential, 2=gaussian)"); + DBL_ARG(a) ; + DBL_ARG(b) ; + OUTPUT(gsl_ran_exppow (r, a, b)); + } + else if (NAME("fdist")) + { + ARGS(2, "nu1, nu2 = degrees of freedom parameters"); + DBL_ARG(nu1) ; + DBL_ARG(nu2) ; + OUTPUT(gsl_ran_fdist (r, nu1, nu2)); + } + else if (NAME("flat")) + { + ARGS(2, "a = lower limit, b = upper limit"); + DBL_ARG(a) ; + DBL_ARG(b) ; + OUTPUT(gsl_ran_flat (r, a, b)); + } + else if (NAME("gamma")) + { + ARGS(2, "a = order, b = scale"); + DBL_ARG(a) ; + DBL_ARG(b) ; + OUTPUT(gsl_ran_gamma (r, a, b)); + } + else if (NAME("gaussian")) + { + ARGS(1, "sigma = standard deviation"); + DBL_ARG(sigma) ; + OUTPUT(gsl_ran_gaussian (r, sigma)); + } + else if (NAME("gaussian-tail")) + { + ARGS(2, "a = lower limit, sigma = standard deviation"); + DBL_ARG(a) ; + DBL_ARG(sigma) ; + OUTPUT(gsl_ran_gaussian_tail (r, a, sigma)); + } + else if (NAME("ugaussian")) + { + ARGS(0, "unit gaussian, no parameters required"); + OUTPUT(gsl_ran_ugaussian (r)); + } + else if (NAME("ugaussian-tail")) + { + ARGS(1, "a = lower limit"); + DBL_ARG(a) ; + OUTPUT(gsl_ran_ugaussian_tail (r, a)); + } + else if (NAME("bivariate-gaussian")) + { + ARGS(3, "sigmax = x std.dev., sigmay = y std.dev., rho = correlation"); + DBL_ARG(sigmax) ; + DBL_ARG(sigmay) ; + DBL_ARG(rho) ; + OUTPUT2(gsl_ran_bivariate_gaussian (r, sigmax, sigmay, rho, &x, &y), + x, y); + } + else if (NAME("dir-2d")) + { + OUTPUT2(gsl_ran_dir_2d (r, &x, &y), x, y); + } + else if (NAME("dir-3d")) + { + OUTPUT3(gsl_ran_dir_3d (r, &x, &y, &z), x, y, z); + } + else if (NAME("dir-nd")) + { + double *xarr; + ARGS(1, "n1 = number of dimensions of hypersphere"); + INT_ARG(n1) ; + xarr = (double *)malloc(n1*sizeof(double)); + + for(i = 0; i < n; i++) { + gsl_ran_dir_nd (r, n1, xarr) ; + for (j = 0; j < n1; j++) { + if (j) putchar(' '); + printf("%g", xarr[j]) ; + } + putchar('\n'); + } ; + + free(xarr); + } + else if (NAME("geometric")) + { + ARGS(1, "p = bernoulli trial probability of success"); + DBL_ARG(p) ; + INT_OUTPUT(gsl_ran_geometric (r, p)); + } + else if (NAME("gumbel1")) + { + ARGS(2, "a = order, b = scale parameter"); + DBL_ARG(a) ; + DBL_ARG(b) ; + OUTPUT(gsl_ran_gumbel1 (r, a, b)); + } + else if (NAME("gumbel2")) + { + ARGS(2, "a = order, b = scale parameter"); + DBL_ARG(a) ; + DBL_ARG(b) ; + OUTPUT(gsl_ran_gumbel2 (r, a, b)); + } + else if (NAME("hypergeometric")) + { + ARGS(3, "n1 = tagged population, n2 = untagged population, t = number of trials"); + INT_ARG(n1) ; + INT_ARG(n2) ; + INT_ARG(t) ; + INT_OUTPUT(gsl_ran_hypergeometric (r, n1, n2, t)); + } + else if (NAME("laplace")) + { + ARGS(1, "a = scale parameter"); + DBL_ARG(a) ; + OUTPUT(gsl_ran_laplace (r, a)); + } + else if (NAME("landau")) + { + ARGS(0, "no arguments required"); + OUTPUT(gsl_ran_landau (r)); + } + else if (NAME("levy")) + { + ARGS(2, "c = scale, a = power (1=cauchy, 2=gaussian)"); + DBL_ARG(c) ; + DBL_ARG(a) ; + OUTPUT(gsl_ran_levy (r, c, a)); + } + else if (NAME("levy-skew")) + { + ARGS(3, "c = scale, a = power (1=cauchy, 2=gaussian), b = skew"); + DBL_ARG(c) ; + DBL_ARG(a) ; + DBL_ARG(b) ; + OUTPUT(gsl_ran_levy_skew (r, c, a, b)); + } + else if (NAME("logarithmic")) + { + ARGS(1, "p = probability"); + DBL_ARG(p) ; + INT_OUTPUT(gsl_ran_logarithmic (r, p)); + } + else if (NAME("logistic")) + { + ARGS(1, "a = scale parameter"); + DBL_ARG(a) ; + OUTPUT(gsl_ran_logistic (r, a)); + } + else if (NAME("lognormal")) + { + ARGS(2, "zeta = location parameter, sigma = scale parameter"); + DBL_ARG(zeta) ; + DBL_ARG(sigma) ; + OUTPUT(gsl_ran_lognormal (r, zeta, sigma)); + } + else if (NAME("negative-binomial")) + { + ARGS(2, "p = probability, a = order"); + DBL_ARG(p) ; + DBL_ARG(a) ; + INT_OUTPUT(gsl_ran_negative_binomial (r, p, a)); + } + else if (NAME("pareto")) + { + ARGS(2, "a = power, b = scale parameter"); + DBL_ARG(a) ; + DBL_ARG(b) ; + OUTPUT(gsl_ran_pareto (r, a, b)); + } + else if (NAME("pascal")) + { + ARGS(2, "p = probability, n = order (integer)"); + DBL_ARG(p) ; + INT_ARG(N) ; + INT_OUTPUT(gsl_ran_pascal (r, p, N)); + } + else if (NAME("poisson")) + { + ARGS(1, "mu = scale parameter"); + DBL_ARG(mu) ; + INT_OUTPUT(gsl_ran_poisson (r, mu)); + } + else if (NAME("rayleigh")) + { + ARGS(1, "sigma = scale parameter"); + DBL_ARG(sigma) ; + OUTPUT(gsl_ran_rayleigh (r, sigma)); + } + else if (NAME("rayleigh-tail")) + { + ARGS(2, "a = lower limit, sigma = scale parameter"); + DBL_ARG(a) ; + DBL_ARG(sigma) ; + OUTPUT(gsl_ran_rayleigh_tail (r, a, sigma)); + } + else if (NAME("tdist")) + { + ARGS(1, "nu = degrees of freedom"); + DBL_ARG(nu) ; + OUTPUT(gsl_ran_tdist (r, nu)); + } + else if (NAME("weibull")) + { + ARGS(2, "a = scale parameter, b = exponent"); + DBL_ARG(a) ; + DBL_ARG(b) ; + OUTPUT(gsl_ran_weibull (r, a, b)); + } + else + { + fprintf(stderr,"Error: unrecognized distribution: %s\n", name) ; + } + + return 0 ; +} + + +void +error (const char * s) +{ + fprintf(stderr, "Error: arguments should be %s\n",s) ; + exit (EXIT_FAILURE) ; +} |