3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
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.
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.
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.
24 #include <gsl/gsl_math.h>
25 #include <gsl/gsl_randist.h>
26 #include <gsl/gsl_rng.h>
27 #include <gsl/gsl_test.h>
28 #include <gsl/gsl_ieee_utils.h>
29 #include <gsl/gsl_integration.h>
33 /* Convient test dimension for multivariant distributions */
37 void testMoments (double (*f) (void), const char *name,
38 double a, double b, double p);
39 void testPDF (double (*f) (void), double (*pdf) (double), const char *name);
40 void testDiscretePDF (double (*f) (void), double (*pdf) (unsigned int),
43 void test_shuffle (void);
44 void test_choose (void);
45 double test_beta (void);
46 double test_beta_pdf (double x);
47 double test_bernoulli (void);
48 double test_bernoulli_pdf (unsigned int n);
50 double test_binomial (void);
51 double test_binomial_pdf (unsigned int n);
52 double test_binomial_large (void);
53 double test_binomial_large_pdf (unsigned int n);
54 double test_binomial_huge (void);
55 double test_binomial_huge_pdf (unsigned int n);
56 double test_binomial0 (void);
57 double test_binomial0_pdf (unsigned int n);
58 double test_binomial1 (void);
59 double test_binomial1_pdf (unsigned int n);
63 double test_binomial_knuth (void);
64 double test_binomial_knuth_pdf (unsigned int n);
65 double test_binomial_large_knuth (void);
66 double test_binomial_large_knuth_pdf (unsigned int n);
67 double test_binomial_huge_knuth (void);
68 double test_binomial_huge_knuth_pdf (unsigned int n);
70 double test_cauchy (void);
71 double test_cauchy_pdf (double x);
72 double test_chisq (void);
73 double test_chisq_pdf (double x);
74 double test_dirichlet (void);
75 double test_dirichlet_pdf (double x);
76 double test_dirichlet_small (void);
77 double test_dirichlet_small_pdf (double x);
78 void test_dirichlet_moments (void);
79 double test_discrete1 (void);
80 double test_discrete1_pdf (unsigned int n);
81 double test_discrete2 (void);
82 double test_discrete2_pdf (unsigned int n);
83 double test_discrete3 (void);
84 double test_discrete3_pdf (unsigned int n);
85 double test_erlang (void);
86 double test_erlang_pdf (double x);
87 double test_exponential (void);
88 double test_exponential_pdf (double x);
89 double test_exppow0 (void);
90 double test_exppow0_pdf (double x);
91 double test_exppow1 (void);
92 double test_exppow1_pdf (double x);
93 double test_exppow1a (void);
94 double test_exppow1a_pdf (double x);
95 double test_exppow2 (void);
96 double test_exppow2_pdf (double x);
97 double test_exppow2a (void);
98 double test_exppow2a_pdf (double x);
99 double test_exppow2b (void);
100 double test_exppow2b_pdf (double x);
101 double test_fdist (void);
102 double test_fdist_pdf (double x);
103 double test_flat (void);
104 double test_flat_pdf (double x);
105 double test_gamma (void);
106 double test_gamma_pdf (double x);
107 double test_gamma1 (void);
108 double test_gamma1_pdf (double x);
109 double test_gamma_int (void);
110 double test_gamma_int_pdf (double x);
111 double test_gamma_large (void);
112 double test_gamma_large_pdf (double x);
113 double test_gamma_small (void);
114 double test_gamma_small_pdf (double x);
115 double test_gamma_mt (void);
116 double test_gamma_mt_pdf (double x);
117 double test_gamma_mt1 (void);
118 double test_gamma_mt1_pdf (double x);
119 double test_gamma_mt_int (void);
120 double test_gamma_mt_int_pdf (double x);
121 double test_gamma_mt_large (void);
122 double test_gamma_mt_large_pdf (double x);
123 double test_gamma_mt_small (void);
124 double test_gamma_mt_small_pdf (double x);
125 double test_gaussian (void);
126 double test_gaussian_pdf (double x);
127 double test_gaussian_ratio_method (void);
128 double test_gaussian_ratio_method_pdf (double x);
129 double test_gaussian_ziggurat (void);
130 double test_gaussian_ziggurat_pdf (double x);
131 double test_gaussian_tail (void);
132 double test_gaussian_tail_pdf (double x);
133 double test_gaussian_tail1 (void);
134 double test_gaussian_tail1_pdf (double x);
135 double test_gaussian_tail2 (void);
136 double test_gaussian_tail2_pdf (double x);
137 double test_ugaussian (void);
138 double test_ugaussian_pdf (double x);
139 double test_ugaussian_ratio_method (void);
140 double test_ugaussian_ratio_method_pdf (double x);
141 double test_ugaussian_tail (void);
142 double test_ugaussian_tail_pdf (double x);
143 double test_bivariate_gaussian1 (void);
144 double test_bivariate_gaussian1_pdf (double x);
145 double test_bivariate_gaussian2 (void);
146 double test_bivariate_gaussian2_pdf (double x);
147 double test_bivariate_gaussian3 (void);
148 double test_bivariate_gaussian3_pdf (double x);
149 double test_bivariate_gaussian4 (void);
150 double test_bivariate_gaussian4_pdf (double x);
151 double test_gumbel1 (void);
152 double test_gumbel1_pdf (double x);
153 double test_gumbel2 (void);
154 double test_gumbel2_pdf (double x);
155 double test_geometric (void);
156 double test_geometric_pdf (unsigned int x);
157 double test_geometric1 (void);
158 double test_geometric1_pdf (unsigned int x);
159 double test_hypergeometric1 (void);
160 double test_hypergeometric1_pdf (unsigned int x);
161 double test_hypergeometric2 (void);
162 double test_hypergeometric2_pdf (unsigned int x);
163 double test_hypergeometric3 (void);
164 double test_hypergeometric3_pdf (unsigned int x);
165 double test_hypergeometric4 (void);
166 double test_hypergeometric4_pdf (unsigned int x);
167 double test_hypergeometric5 (void);
168 double test_hypergeometric5_pdf (unsigned int x);
169 double test_hypergeometric6 (void);
170 double test_hypergeometric6_pdf (unsigned int x);
171 double test_landau (void);
172 double test_landau_pdf (double x);
173 double test_levy1 (void);
174 double test_levy1_pdf (double x);
175 double test_levy2 (void);
176 double test_levy2_pdf (double x);
177 double test_levy1a (void);
178 double test_levy1a_pdf (double x);
179 double test_levy2a (void);
180 double test_levy2a_pdf (double x);
181 double test_levy_skew1 (void);
182 double test_levy_skew1_pdf (double x);
183 double test_levy_skew2 (void);
184 double test_levy_skew2_pdf (double x);
185 double test_levy_skew1a (void);
186 double test_levy_skew1a_pdf (double x);
187 double test_levy_skew2a (void);
188 double test_levy_skew2a_pdf (double x);
189 double test_levy_skew1b (void);
190 double test_levy_skew1b_pdf (double x);
191 double test_levy_skew2b (void);
192 double test_levy_skew2b_pdf (double x);
193 double test_logistic (void);
194 double test_logistic_pdf (double x);
195 double test_lognormal (void);
196 double test_lognormal_pdf (double x);
197 double test_logarithmic (void);
198 double test_logarithmic_pdf (unsigned int n);
199 double test_multinomial (void);
200 double test_multinomial_pdf (unsigned int n);
201 double test_multinomial_large (void);
202 double test_multinomial_large_pdf (unsigned int n);
203 void test_multinomial_moments (void);
204 double test_negative_binomial (void);
205 double test_negative_binomial_pdf (unsigned int n);
206 double test_pascal (void);
207 double test_pascal_pdf (unsigned int n);
208 double test_pareto (void);
209 double test_pareto_pdf (double x);
210 double test_poisson (void);
211 double test_poisson_pdf (unsigned int x);
212 double test_poisson_large (void);
213 double test_poisson_large_pdf (unsigned int x);
214 double test_dir2d (void);
215 double test_dir2d_pdf (double x);
216 double test_dir2d_trig_method (void);
217 double test_dir2d_trig_method_pdf (double x);
218 double test_dir3dxy (void);
219 double test_dir3dxy_pdf (double x);
220 double test_dir3dyz (void);
221 double test_dir3dyz_pdf (double x);
222 double test_dir3dzx (void);
223 double test_dir3dzx_pdf (double x);
224 double test_rayleigh (void);
225 double test_rayleigh_pdf (double x);
226 double test_rayleigh_tail (void);
227 double test_rayleigh_tail_pdf (double x);
228 double test_tdist1 (void);
229 double test_tdist1_pdf (double x);
230 double test_tdist2 (void);
231 double test_tdist2_pdf (double x);
232 double test_laplace (void);
233 double test_laplace_pdf (double x);
234 double test_weibull (void);
235 double test_weibull_pdf (double x);
236 double test_weibull1 (void);
237 double test_weibull1_pdf (double x);
241 static gsl_ran_discrete_t *g1 = NULL;
242 static gsl_ran_discrete_t *g2 = NULL;
243 static gsl_ran_discrete_t *g3 = NULL;
248 gsl_ieee_env_setup ();
250 gsl_rng_env_setup ();
251 r_global = gsl_rng_alloc (gsl_rng_default);
253 #define FUNC(x) test_ ## x, "test gsl_ran_" #x
254 #define FUNC2(x) test_ ## x, test_ ## x ## _pdf, "test gsl_ran_" #x
259 testMoments (FUNC (ugaussian), 0.0, 100.0, 0.5);
260 testMoments (FUNC (ugaussian), -1.0, 1.0, 0.6826895);
261 testMoments (FUNC (ugaussian), 3.0, 3.5, 0.0011172689);
262 testMoments (FUNC (ugaussian_tail), 3.0, 3.5, 0.0011172689 / 0.0013498981);
263 testMoments (FUNC (exponential), 0.0, 1.0, 1 - exp (-0.5));
264 testMoments (FUNC (cauchy), 0.0, 10000.0, 0.5);
266 testMoments (FUNC (discrete1), -0.5, 0.5, 0.59);
267 testMoments (FUNC (discrete1), 0.5, 1.5, 0.40);
268 testMoments (FUNC (discrete1), 1.5, 3.5, 0.01);
270 testMoments (FUNC (discrete2), -0.5, 0.5, 1.0/45.0 );
271 testMoments (FUNC (discrete2), 8.5, 9.5, 0 );
273 testMoments (FUNC (discrete3), -0.5, 0.5, 0.05 );
274 testMoments (FUNC (discrete3), 0.5, 1.5, 0.05 );
275 testMoments (FUNC (discrete3), -0.5, 9.5, 0.5 );
277 test_dirichlet_moments ();
278 test_multinomial_moments ();
280 testPDF (FUNC2 (beta));
281 testPDF (FUNC2 (cauchy));
282 testPDF (FUNC2 (chisq));
283 testPDF (FUNC2 (dirichlet));
284 testPDF (FUNC2 (dirichlet_small));
285 testPDF (FUNC2 (erlang));
286 testPDF (FUNC2 (exponential));
288 testPDF (FUNC2 (exppow0));
289 testPDF (FUNC2 (exppow1));
290 testPDF (FUNC2 (exppow1a));
291 testPDF (FUNC2 (exppow2));
292 testPDF (FUNC2 (exppow2a));
293 testPDF (FUNC2 (exppow2b));
295 testPDF (FUNC2 (fdist));
296 testPDF (FUNC2 (flat));
297 testPDF (FUNC2 (gamma));
298 testPDF (FUNC2 (gamma1));
299 testPDF (FUNC2 (gamma_int));
300 testPDF (FUNC2 (gamma_large));
301 testPDF (FUNC2 (gamma_small));
302 testPDF (FUNC2 (gamma_mt));
303 testPDF (FUNC2 (gamma_mt1));
304 testPDF (FUNC2 (gamma_mt_int));
305 testPDF (FUNC2 (gamma_mt_large));
306 testPDF (FUNC2 (gamma_mt_small));
307 testPDF (FUNC2 (gaussian));
308 testPDF (FUNC2 (gaussian_ratio_method));
309 testPDF (FUNC2 (gaussian_ziggurat));
310 testPDF (FUNC2 (ugaussian));
311 testPDF (FUNC2 (ugaussian_ratio_method));
312 testPDF (FUNC2 (gaussian_tail));
313 testPDF (FUNC2 (gaussian_tail1));
314 testPDF (FUNC2 (gaussian_tail2));
315 testPDF (FUNC2 (ugaussian_tail));
317 testPDF (FUNC2 (bivariate_gaussian1));
318 testPDF (FUNC2 (bivariate_gaussian2));
319 testPDF (FUNC2 (bivariate_gaussian3));
320 testPDF (FUNC2 (bivariate_gaussian4));
322 testPDF (FUNC2 (gumbel1));
323 testPDF (FUNC2 (gumbel2));
324 testPDF (FUNC2 (landau));
325 testPDF (FUNC2 (levy1));
326 testPDF (FUNC2 (levy2));
327 testPDF (FUNC2 (levy1a));
328 testPDF (FUNC2 (levy2a));
329 testPDF (FUNC2 (levy_skew1));
330 testPDF (FUNC2 (levy_skew2));
331 testPDF (FUNC2 (levy_skew1a));
332 testPDF (FUNC2 (levy_skew2a));
333 testPDF (FUNC2 (levy_skew1b));
334 testPDF (FUNC2 (levy_skew2b));
335 testPDF (FUNC2 (logistic));
336 testPDF (FUNC2 (lognormal));
337 testPDF (FUNC2 (pareto));
338 testPDF (FUNC2 (rayleigh));
339 testPDF (FUNC2 (rayleigh_tail));
340 testPDF (FUNC2 (tdist1));
341 testPDF (FUNC2 (tdist2));
342 testPDF (FUNC2 (laplace));
343 testPDF (FUNC2 (weibull));
344 testPDF (FUNC2 (weibull1));
346 testPDF (FUNC2 (dir2d));
347 testPDF (FUNC2 (dir2d_trig_method));
348 testPDF (FUNC2 (dir3dxy));
349 testPDF (FUNC2 (dir3dyz));
350 testPDF (FUNC2 (dir3dzx));
352 testDiscretePDF (FUNC2 (discrete1));
353 testDiscretePDF (FUNC2 (discrete2));
354 testDiscretePDF (FUNC2 (discrete3));
355 testDiscretePDF (FUNC2 (poisson));
356 testDiscretePDF (FUNC2 (poisson_large));
357 testDiscretePDF (FUNC2 (bernoulli));
358 testDiscretePDF (FUNC2 (binomial));
359 testDiscretePDF (FUNC2 (binomial0));
360 testDiscretePDF (FUNC2 (binomial1));
361 testDiscretePDF (FUNC2 (binomial_knuth));
362 testDiscretePDF (FUNC2 (binomial_large));
363 testDiscretePDF (FUNC2 (binomial_large_knuth));
364 testDiscretePDF (FUNC2 (binomial_huge));
365 testDiscretePDF (FUNC2 (binomial_huge_knuth));
366 testDiscretePDF (FUNC2 (geometric));
367 testDiscretePDF (FUNC2 (geometric1));
368 testDiscretePDF (FUNC2 (hypergeometric1));
369 testDiscretePDF (FUNC2 (hypergeometric2));
370 testDiscretePDF (FUNC2 (hypergeometric3));
371 testDiscretePDF (FUNC2 (hypergeometric4));
372 testDiscretePDF (FUNC2 (hypergeometric5));
373 testDiscretePDF (FUNC2 (hypergeometric6));
374 testDiscretePDF (FUNC2 (logarithmic));
375 testDiscretePDF (FUNC2 (multinomial));
376 testDiscretePDF (FUNC2 (multinomial_large));
377 testDiscretePDF (FUNC2 (negative_binomial));
378 testDiscretePDF (FUNC2 (pascal));
380 gsl_rng_free (r_global);
381 gsl_ran_discrete_free (g1);
382 gsl_ran_discrete_free (g2);
383 gsl_ran_discrete_free (g3);
385 exit (gsl_test_summary ());
391 double count[10][10];
392 int x[10] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 };
393 int i, j, status = 0;
395 for (i = 0; i < 10; i++)
397 for (j = 0; j < 10; j++)
403 for (i = 0; i < N; i++)
405 for (j = 0; j < 10; j++)
408 gsl_ran_shuffle (r_global, x, 10, sizeof (int));
410 for (j = 0; j < 10; j++)
414 for (i = 0; i < 10; i++)
416 for (j = 0; j < 10; j++)
418 double expected = N / 10.0;
419 double d = fabs (count[i][j] - expected);
420 double sigma = d / sqrt (expected);
421 if (sigma > 5 && d > 1)
425 "gsl_ran_shuffle %d,%d (%g observed vs %g expected)",
426 i, j, count[i][j] / N, 0.1);
431 gsl_test (status, "gsl_ran_shuffle on {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}");
439 int x[10] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 };
440 int y[3] = { 0, 1, 2 };
441 int i, j, status = 0;
443 for (i = 0; i < 10; i++)
448 for (i = 0; i < N; i++)
450 for (j = 0; j < 10; j++)
453 gsl_ran_choose (r_global, y, 3, x, 10, sizeof (int));
455 for (j = 0; j < 3; j++)
459 for (i = 0; i < 10; i++)
461 double expected = 3.0 * N / 10.0;
462 double d = fabs (count[i] - expected);
463 double sigma = d / sqrt (expected);
464 if (sigma > 5 && d > 1)
468 "gsl_ran_choose %d (%g observed vs %g expected)",
469 i, count[i] / N, 0.1);
473 gsl_test (status, "gsl_ran_choose (3) on {0, 1, 2, 3, 4, 5, 6, 7, 8, 9}");
481 testMoments (double (*f) (void), const char *name,
482 double a, double b, double p)
485 double count = 0, expected, sigma;
488 for (i = 0; i < N; i++)
496 sigma = (expected > 0) ? fabs (count - expected) / sqrt (expected) : fabs(count - expected);
498 status = (sigma > 3);
500 gsl_test (status, "%s [%g,%g] (%g observed vs %g expected)",
501 name, a, b, count / N, p);
506 typedef double pdf_func(double);
509 wrapper_function (double x, void *params)
511 pdf_func * pdf = (pdf_func *)params;
516 integrate (pdf_func * pdf, double a, double b)
518 double result, abserr;
521 gsl_integration_workspace * w = gsl_integration_workspace_alloc (n);
522 f.function = &wrapper_function;
523 f.params = (void *)pdf;
524 gsl_integration_qags (&f, a, b, 1e-16, 1e-4, n, w, &result, &abserr);
525 gsl_integration_workspace_free (w);
531 testPDF (double (*f) (void), double (*pdf) (double), const char *name)
533 double count[BINS], edge[BINS], p[BINS];
534 double a = -5.0, b = +5.0;
535 double dx = (b - a) / BINS;
537 double total = 0, mean;
538 int i, j, status = 0, status_i = 0;
540 for (i = 0; i < BINS; i++)
546 for (i = 0; i < N; i++)
554 double u = (r - a) / dx;
555 double f = modf(u, &bin);
565 /* Sort out where the hits on the edges should go */
567 for (i = 0; i < BINS; i++)
569 /* If the bin above is empty, its lower edge hits belong in the
572 if (i + 1 < BINS && count[i+1] == 0) {
573 count[i] += edge[i+1];
582 gsl_test (!gsl_finite(mean), "%s, finite mean, observed %g", name, mean);
584 for (i = 0; i < BINS; i++)
586 /* Compute an approximation to the integral of p(x) from x to
587 x+dx using Simpson's rule */
589 double x = a + i * dx;
591 if (fabs (x) < 1e-10) /* hit the origin exactly */
594 p[i] = integrate (pdf, x, x+dx);
597 for (i = 0; i < BINS; i++)
599 double x = a + i * dx;
600 double d = fabs (count[i] - N * p[i]);
603 double s = d / sqrt (N * p[i]);
604 status_i = (s > 5) && (d > 2);
608 status_i = (count[i] != 0);
612 gsl_test (status_i, "%s [%g,%g) (%g/%d=%g observed vs %g expected)",
613 name, x, x + dx, count[i], N, count[i] / N, p[i]);
617 gsl_test (status, "%s, sampling against pdf over range [%g,%g) ",
622 testDiscretePDF (double (*f) (void), double (*pdf) (unsigned int),
625 double count[BINS], p[BINS];
627 int status = 0, status_i = 0;
629 for (i = 0; i < BINS; i++)
632 for (i = 0; i < N; i++)
634 int r = (int) (f ());
635 if (r >= 0 && r < BINS)
639 for (i = 0; i < BINS; i++)
642 for (i = 0; i < BINS; i++)
644 double d = fabs (count[i] - N * p[i]);
647 double s = d / sqrt (N * p[i]);
648 status_i = (s > 5) && (d > 1);
652 status_i = (count[i] != 0);
656 gsl_test (status_i, "%s i=%d (%g observed vs %g expected)",
657 name, i, count[i] / N, p[i]);
661 gsl_test (status, "%s, sampling against pdf over range [%d,%d) ",
670 return gsl_ran_beta (r_global, 2.0, 3.0);
674 test_beta_pdf (double x)
676 return gsl_ran_beta_pdf (x, 2.0, 3.0);
680 test_bernoulli (void)
682 return gsl_ran_bernoulli (r_global, 0.3);
686 test_bernoulli_pdf (unsigned int n)
688 return gsl_ran_bernoulli_pdf (n, 0.3);
694 return gsl_ran_binomial (r_global, 0.3, 5);
698 test_binomial_pdf (unsigned int n)
700 return gsl_ran_binomial_pdf (n, 0.3, 5);
704 test_binomial0 (void)
706 return gsl_ran_binomial (r_global, 0, 8);
710 test_binomial0_pdf (unsigned int n)
712 return gsl_ran_binomial_pdf (n, 0, 8);
716 test_binomial1 (void)
718 return gsl_ran_binomial (r_global, 1, 8);
722 test_binomial1_pdf (unsigned int n)
724 return gsl_ran_binomial_pdf (n, 1, 8);
728 test_binomial_knuth (void)
730 return gsl_ran_binomial_knuth (r_global, 0.3, 5);
734 test_binomial_knuth_pdf (unsigned int n)
736 return gsl_ran_binomial_pdf (n, 0.3, 5);
741 test_binomial_large (void)
743 return gsl_ran_binomial (r_global, 0.3, 55);
747 test_binomial_large_pdf (unsigned int n)
749 return gsl_ran_binomial_pdf (n, 0.3, 55);
753 test_binomial_large_knuth (void)
755 return gsl_ran_binomial_knuth (r_global, 0.3, 55);
759 test_binomial_large_knuth_pdf (unsigned int n)
761 return gsl_ran_binomial_pdf (n, 0.3, 55);
766 test_binomial_huge (void)
768 return gsl_ran_binomial (r_global, 0.3, 5500);
772 test_binomial_huge_pdf (unsigned int n)
774 return gsl_ran_binomial_pdf (n, 0.3, 5500);
778 test_binomial_huge_knuth (void)
780 return gsl_ran_binomial_knuth (r_global, 0.3, 5500);
784 test_binomial_huge_knuth_pdf (unsigned int n)
786 return gsl_ran_binomial_pdf (n, 0.3, 5500);
792 return gsl_ran_cauchy (r_global, 2.0);
796 test_cauchy_pdf (double x)
798 return gsl_ran_cauchy_pdf (x, 2.0);
804 return gsl_ran_chisq (r_global, 13.0);
808 test_chisq_pdf (double x)
810 return gsl_ran_chisq_pdf (x, 13.0);
816 double x = 0, y = 0, theta;
817 gsl_ran_dir_2d (r_global, &x, &y);
818 theta = atan2 (x, y);
823 test_dir2d_pdf (double x)
825 if (x > -M_PI && x <= M_PI)
827 return 1 / (2 * M_PI);
836 test_dir2d_trig_method (void)
838 double x = 0, y = 0, theta;
839 gsl_ran_dir_2d_trig_method (r_global, &x, &y);
840 theta = atan2 (x, y);
845 test_dir2d_trig_method_pdf (double x)
847 if (x > -M_PI && x <= M_PI)
849 return 1 / (2 * M_PI);
860 double x = 0, y = 0, z = 0, theta;
861 gsl_ran_dir_3d (r_global, &x, &y, &z);
862 theta = atan2 (x, y);
867 test_dir3dxy_pdf (double x)
869 if (x > -M_PI && x <= M_PI)
871 return 1 / (2 * M_PI);
882 double x = 0, y = 0, z = 0, theta;
883 gsl_ran_dir_3d (r_global, &x, &y, &z);
884 theta = atan2 (y, z);
889 test_dir3dyz_pdf (double x)
891 if (x > -M_PI && x <= M_PI)
893 return 1 / (2 * M_PI);
904 double x = 0, y = 0, z = 0, theta;
905 gsl_ran_dir_3d (r_global, &x, &y, &z);
906 theta = atan2 (z, x);
911 test_dir3dzx_pdf (double x)
913 if (x > -M_PI && x <= M_PI)
915 return 1 / (2 * M_PI);
924 test_dirichlet (void)
926 /* This is a bit of a lame test, since when K=2, the Dirichlet distribution
927 becomes a beta distribution */
929 double alpha[2] = { 2.5, 5.0 };
930 double theta[2] = { 0.0, 0.0 };
932 gsl_ran_dirichlet (r_global, K, alpha, theta);
938 test_dirichlet_pdf (double x)
941 double alpha[2] = { 2.5, 5.0 };
944 if (x <= 0.0 || x >= 1.0)
945 return 0.0; /* Out of range */
950 return gsl_ran_dirichlet_pdf (K, alpha, theta);
955 test_dirichlet_small (void)
958 double alpha[2] = { 2.5e-3, 5.0e-3};
959 double theta[2] = { 0.0, 0.0 };
961 gsl_ran_dirichlet (r_global, K, alpha, theta);
967 test_dirichlet_small_pdf (double x)
970 double alpha[2] = { 2.5e-3, 5.0e-3 };
973 if (x <= 0.0 || x >= 1.0)
974 return 0.0; /* Out of range */
979 return gsl_ran_dirichlet_pdf (K, alpha, theta);
983 /* Check that the observed means of the Dirichlet variables are
984 within reasonable statistical errors of their correct values. */
986 #define DIRICHLET_K 10
989 test_dirichlet_moments (void)
991 double alpha[DIRICHLET_K];
992 double theta[DIRICHLET_K];
993 double theta_sum[DIRICHLET_K];
995 double alpha_sum = 0.0;
996 double mean, obs_mean, sd, sigma;
999 for (k = 0; k < DIRICHLET_K; k++)
1001 alpha[k] = gsl_ran_exponential (r_global, 0.1);
1002 alpha_sum += alpha[k];
1006 for (n = 0; n < N; n++)
1008 gsl_ran_dirichlet (r_global, DIRICHLET_K, alpha, theta);
1009 for (k = 0; k < DIRICHLET_K; k++)
1010 theta_sum[k] += theta[k];
1013 for (k = 0; k < DIRICHLET_K; k++)
1015 mean = alpha[k] / alpha_sum;
1017 sqrt ((alpha[k] * (1. - alpha[k] / alpha_sum)) /
1018 (alpha_sum * (alpha_sum + 1.)));
1019 obs_mean = theta_sum[k] / N;
1020 sigma = sqrt ((double) N) * fabs (mean - obs_mean) / sd;
1022 status = (sigma > 3.0);
1025 "test gsl_ran_dirichlet: mean (%g observed vs %g expected)",
1031 /* Check that the observed means of the multinomial variables are
1032 within reasonable statistical errors of their correct values. */
1035 test_multinomial_moments (void)
1037 const unsigned int sum_n = 100;
1039 const double p[MULTI_DIM] ={ 0.2, 0.20, 0.17, 0.14, 0.12,
1040 0.07, 0.05, 0.02, 0.02, 0.01 };
1042 unsigned int x[MULTI_DIM];
1043 double x_sum[MULTI_DIM];
1045 double mean, obs_mean, sd, sigma;
1048 for (k = 0; k < MULTI_DIM; k++)
1051 for (n = 0; n < N; n++)
1053 gsl_ran_multinomial (r_global, MULTI_DIM, sum_n, p, x);
1054 for (k = 0; k < MULTI_DIM; k++)
1058 for (k = 0; k < MULTI_DIM; k++)
1060 mean = p[k] * sum_n;
1061 sd = p[k] * (1.-p[k]) * sum_n;
1063 obs_mean = x_sum[k] / N;
1064 sigma = sqrt ((double) N) * fabs (mean - obs_mean) / sd;
1066 status = (sigma > 3.0);
1069 "test gsl_ran_multinomial: mean (%g observed vs %g expected)",
1076 test_discrete1 (void)
1078 static double P[3] = { 0.59, 0.4, 0.01 };
1081 g1 = gsl_ran_discrete_preproc (3, P);
1083 return gsl_ran_discrete (r_global, g1);
1087 test_discrete1_pdf (unsigned int n)
1089 return gsl_ran_discrete_pdf ((size_t) n, g1);
1093 test_discrete2 (void)
1095 static double P[10] = { 1, 9, 3, 4, 5, 8, 6, 7, 2, 0 };
1098 g2 = gsl_ran_discrete_preproc (10, P);
1100 return gsl_ran_discrete (r_global, g2);
1104 test_discrete2_pdf (unsigned int n)
1106 return gsl_ran_discrete_pdf ((size_t) n, g2);
1109 test_discrete3 (void)
1111 static double P[20];
1114 for (i=0; i<20; ++i) P[i]=1.0/20;
1115 g3 = gsl_ran_discrete_preproc (20, P);
1117 return gsl_ran_discrete (r_global, g3);
1121 test_discrete3_pdf (unsigned int n)
1123 return gsl_ran_discrete_pdf ((size_t) n, g3);
1130 return gsl_ran_erlang (r_global, 3.0, 4.0);
1134 test_erlang_pdf (double x)
1136 return gsl_ran_erlang_pdf (x, 3.0, 4.0);
1140 test_exponential (void)
1142 return gsl_ran_exponential (r_global, 2.0);
1146 test_exponential_pdf (double x)
1148 return gsl_ran_exponential_pdf (x, 2.0);
1154 return gsl_ran_exppow (r_global, 3.7, 0.3);
1158 test_exppow0_pdf (double x)
1160 return gsl_ran_exppow_pdf (x, 3.7, 0.3);
1166 return gsl_ran_exppow (r_global, 3.7, 1.0);
1170 test_exppow1_pdf (double x)
1172 return gsl_ran_exppow_pdf (x, 3.7, 1.0);
1176 test_exppow1a (void)
1178 return gsl_ran_exppow (r_global, 3.7, 1.9);
1182 test_exppow1a_pdf (double x)
1184 return gsl_ran_exppow_pdf (x, 3.7, 1.9);
1190 return gsl_ran_exppow (r_global, 3.7, 2.0);
1194 test_exppow2_pdf (double x)
1196 return gsl_ran_exppow_pdf (x, 3.7, 2.0);
1201 test_exppow2a (void)
1203 return gsl_ran_exppow (r_global, 3.7, 3.5);
1207 test_exppow2a_pdf (double x)
1209 return gsl_ran_exppow_pdf (x, 3.7, 3.5);
1213 test_exppow2b (void)
1215 return gsl_ran_exppow (r_global, 3.7, 7.5);
1219 test_exppow2b_pdf (double x)
1221 return gsl_ran_exppow_pdf (x, 3.7, 7.5);
1227 return gsl_ran_fdist (r_global, 3.0, 4.0);
1231 test_fdist_pdf (double x)
1233 return gsl_ran_fdist_pdf (x, 3.0, 4.0);
1239 return gsl_ran_flat (r_global, 3.0, 4.0);
1243 test_flat_pdf (double x)
1245 return gsl_ran_flat_pdf (x, 3.0, 4.0);
1251 return gsl_ran_gamma (r_global, 2.5, 2.17);
1255 test_gamma_pdf (double x)
1257 return gsl_ran_gamma_pdf (x, 2.5, 2.17);
1263 return gsl_ran_gamma (r_global, 1.0, 2.17);
1267 test_gamma1_pdf (double x)
1269 return gsl_ran_gamma_pdf (x, 1.0, 2.17);
1274 test_gamma_int (void)
1276 return gsl_ran_gamma (r_global, 10.0, 2.17);
1280 test_gamma_int_pdf (double x)
1282 return gsl_ran_gamma_pdf (x, 10.0, 2.17);
1287 test_gamma_large (void)
1289 return gsl_ran_gamma (r_global, 20.0, 2.17);
1293 test_gamma_large_pdf (double x)
1295 return gsl_ran_gamma_pdf (x, 20.0, 2.17);
1299 test_gamma_small (void)
1301 return gsl_ran_gamma (r_global, 0.92, 2.17);
1305 test_gamma_small_pdf (double x)
1307 return gsl_ran_gamma_pdf (x, 0.92, 2.17);
1312 test_gamma_mt (void)
1314 return gsl_ran_gamma_mt (r_global, 2.5, 2.17);
1318 test_gamma_mt_pdf (double x)
1320 return gsl_ran_gamma_pdf (x, 2.5, 2.17);
1324 test_gamma_mt1 (void)
1326 return gsl_ran_gamma_mt (r_global, 1.0, 2.17);
1330 test_gamma_mt1_pdf (double x)
1332 return gsl_ran_gamma_pdf (x, 1.0, 2.17);
1337 test_gamma_mt_int (void)
1339 return gsl_ran_gamma_mt (r_global, 10.0, 2.17);
1343 test_gamma_mt_int_pdf (double x)
1345 return gsl_ran_gamma_pdf (x, 10.0, 2.17);
1350 test_gamma_mt_large (void)
1352 return gsl_ran_gamma_mt (r_global, 20.0, 2.17);
1356 test_gamma_mt_large_pdf (double x)
1358 return gsl_ran_gamma_pdf (x, 20.0, 2.17);
1363 test_gamma_mt_small (void)
1365 return gsl_ran_gamma_mt (r_global, 0.92, 2.17);
1369 test_gamma_mt_small_pdf (double x)
1371 return gsl_ran_gamma_pdf (x, 0.92, 2.17);
1376 test_gaussian (void)
1378 return gsl_ran_gaussian (r_global, 3.0);
1382 test_gaussian_pdf (double x)
1384 return gsl_ran_gaussian_pdf (x, 3.0);
1388 test_gaussian_ratio_method (void)
1390 return gsl_ran_gaussian_ratio_method (r_global, 3.0);
1394 test_gaussian_ratio_method_pdf (double x)
1396 return gsl_ran_gaussian_pdf (x, 3.0);
1400 test_gaussian_ziggurat (void)
1402 return gsl_ran_gaussian_ziggurat (r_global, 3.12);
1406 test_gaussian_ziggurat_pdf (double x)
1408 return gsl_ran_gaussian_pdf (x, 3.12);
1412 test_gaussian_tail (void)
1414 return gsl_ran_gaussian_tail (r_global, 1.7, 0.25);
1418 test_gaussian_tail_pdf (double x)
1420 return gsl_ran_gaussian_tail_pdf (x, 1.7, 0.25);
1424 test_gaussian_tail1 (void)
1426 return gsl_ran_gaussian_tail (r_global, -1.7, 5.0);
1430 test_gaussian_tail1_pdf (double x)
1432 return gsl_ran_gaussian_tail_pdf (x, -1.7, 5.0);
1436 test_gaussian_tail2 (void)
1438 return gsl_ran_gaussian_tail (r_global, 0.1, 2.0);
1442 test_gaussian_tail2_pdf (double x)
1444 return gsl_ran_gaussian_tail_pdf (x, 0.1, 2.0);
1449 test_ugaussian (void)
1451 return gsl_ran_ugaussian (r_global);
1455 test_ugaussian_pdf (double x)
1457 return gsl_ran_ugaussian_pdf (x);
1461 test_ugaussian_ratio_method (void)
1463 return gsl_ran_ugaussian_ratio_method (r_global);
1467 test_ugaussian_ratio_method_pdf (double x)
1469 return gsl_ran_ugaussian_pdf (x);
1473 test_ugaussian_tail (void)
1475 return gsl_ran_ugaussian_tail (r_global, 3.0);
1479 test_ugaussian_tail_pdf (double x)
1481 return gsl_ran_ugaussian_tail_pdf (x, 3.0);
1485 test_bivariate_gaussian1 (void)
1487 double x = 0, y = 0;
1488 gsl_ran_bivariate_gaussian (r_global, 3.0, 2.0, 0.3, &x, &y);
1493 test_bivariate_gaussian1_pdf (double x)
1495 return gsl_ran_gaussian_pdf (x, 3.0);
1499 test_bivariate_gaussian2 (void)
1501 double x = 0, y = 0;
1502 gsl_ran_bivariate_gaussian (r_global, 3.0, 2.0, 0.3, &x, &y);
1507 test_bivariate_gaussian2_pdf (double y)
1511 double a = -10, b = 10, dx = (b - a) / n;
1512 for (i = 0; i < n; i++)
1514 double x = a + i * dx;
1515 sum += gsl_ran_bivariate_gaussian_pdf (x, y, 3.0, 2.0, 0.3) * dx;
1522 test_bivariate_gaussian3 (void)
1524 double x = 0, y = 0;
1525 gsl_ran_bivariate_gaussian (r_global, 3.0, 2.0, 0.3, &x, &y);
1530 test_bivariate_gaussian3_pdf (double x)
1532 double sx = 3.0, sy = 2.0, r = 0.3;
1533 double su = (sx + r * sy);
1534 double sv = sy * sqrt (1 - r * r);
1535 double sigma = sqrt (su * su + sv * sv);
1537 return gsl_ran_gaussian_pdf (x, sigma);
1541 test_bivariate_gaussian4 (void)
1543 double x = 0, y = 0;
1544 gsl_ran_bivariate_gaussian (r_global, 3.0, 2.0, -0.5, &x, &y);
1549 test_bivariate_gaussian4_pdf (double x)
1551 double sx = 3.0, sy = 2.0, r = -0.5;
1552 double su = (sx + r * sy);
1553 double sv = sy * sqrt (1 - r * r);
1554 double sigma = sqrt (su * su + sv * sv);
1556 return gsl_ran_gaussian_pdf (x, sigma);
1561 test_geometric (void)
1563 return gsl_ran_geometric (r_global, 0.5);
1567 test_geometric_pdf (unsigned int n)
1569 return gsl_ran_geometric_pdf (n, 0.5);
1573 test_geometric1 (void)
1575 return gsl_ran_geometric (r_global, 1.0);
1579 test_geometric1_pdf (unsigned int n)
1581 return gsl_ran_geometric_pdf (n, 1.0);
1585 test_hypergeometric1 (void)
1587 return gsl_ran_hypergeometric (r_global, 5, 7, 4);
1591 test_hypergeometric1_pdf (unsigned int n)
1593 return gsl_ran_hypergeometric_pdf (n, 5, 7, 4);
1598 test_hypergeometric2 (void)
1600 return gsl_ran_hypergeometric (r_global, 5, 7, 11);
1604 test_hypergeometric2_pdf (unsigned int n)
1606 return gsl_ran_hypergeometric_pdf (n, 5, 7, 11);
1610 test_hypergeometric3 (void)
1612 return gsl_ran_hypergeometric (r_global, 5, 7, 1);
1616 test_hypergeometric3_pdf (unsigned int n)
1618 return gsl_ran_hypergeometric_pdf (n, 5, 7, 1);
1622 test_hypergeometric4 (void)
1624 return gsl_ran_hypergeometric (r_global, 5, 7, 20);
1628 test_hypergeometric4_pdf (unsigned int n)
1630 return gsl_ran_hypergeometric_pdf (n, 5, 7, 20);
1634 test_hypergeometric5 (void)
1636 return gsl_ran_hypergeometric (r_global, 2, 7, 5);
1640 test_hypergeometric5_pdf (unsigned int n)
1642 return gsl_ran_hypergeometric_pdf (n, 2, 7, 5);
1647 test_hypergeometric6 (void)
1649 return gsl_ran_hypergeometric (r_global, 2, 10, 3);
1653 test_hypergeometric6_pdf (unsigned int n)
1655 return gsl_ran_hypergeometric_pdf (n, 2, 10, 3);
1664 return gsl_ran_gumbel1 (r_global, 3.12, 4.56);
1668 test_gumbel1_pdf (double x)
1670 return gsl_ran_gumbel1_pdf (x, 3.12, 4.56);
1676 return gsl_ran_gumbel2 (r_global, 3.12, 4.56);
1680 test_gumbel2_pdf (double x)
1682 return gsl_ran_gumbel2_pdf (x, 3.12, 4.56);
1688 return gsl_ran_landau (r_global);
1692 test_landau_pdf (double x)
1694 return gsl_ran_landau_pdf (x);
1700 return gsl_ran_levy (r_global, 5.0, 1.0);
1704 test_levy1_pdf (double x)
1706 return gsl_ran_cauchy_pdf (x, 5.0);
1712 return gsl_ran_levy (r_global, 5.0, 2.0);
1716 test_levy2_pdf (double x)
1718 return gsl_ran_gaussian_pdf (x, sqrt (2.0) * 5.0);
1724 return gsl_ran_levy (r_global, 5.0, 1.01);
1728 test_levy1a_pdf (double x)
1730 return gsl_ran_cauchy_pdf (x, 5.0);
1736 return gsl_ran_levy (r_global, 5.0, 1.99);
1740 test_levy2a_pdf (double x)
1742 return gsl_ran_gaussian_pdf (x, sqrt (2.0) * 5.0);
1747 test_levy_skew1 (void)
1749 return gsl_ran_levy_skew (r_global, 5.0, 1.0, 0.0);
1753 test_levy_skew1_pdf (double x)
1755 return gsl_ran_cauchy_pdf (x, 5.0);
1759 test_levy_skew2 (void)
1761 return gsl_ran_levy_skew (r_global, 5.0, 2.0, 0.0);
1765 test_levy_skew2_pdf (double x)
1767 return gsl_ran_gaussian_pdf (x, sqrt (2.0) * 5.0);
1771 test_levy_skew1a (void)
1773 return gsl_ran_levy_skew (r_global, 5.0, 1.01, 0.0);
1777 test_levy_skew1a_pdf (double x)
1779 return gsl_ran_cauchy_pdf (x, 5.0);
1783 test_levy_skew2a (void)
1785 return gsl_ran_levy_skew (r_global, 5.0, 1.99, 0.0);
1789 test_levy_skew2a_pdf (double x)
1791 return gsl_ran_gaussian_pdf (x, sqrt (2.0) * 5.0);
1795 test_levy_skew1b (void)
1797 return gsl_ran_levy_skew (r_global, 5.0, 1.01, 0.001);
1801 test_levy_skew1b_pdf (double x)
1803 return gsl_ran_cauchy_pdf (x, 5.0);
1807 test_levy_skew2b (void)
1809 return gsl_ran_levy_skew (r_global, 5.0, 1.99, 0.001);
1813 test_levy_skew2b_pdf (double x)
1815 return gsl_ran_gaussian_pdf (x, sqrt (2.0) * 5.0);
1820 test_logistic (void)
1822 return gsl_ran_logistic (r_global, 3.1);
1826 test_logistic_pdf (double x)
1828 return gsl_ran_logistic_pdf (x, 3.1);
1832 test_logarithmic (void)
1834 return gsl_ran_logarithmic (r_global, 0.4);
1838 test_logarithmic_pdf (unsigned int n)
1840 return gsl_ran_logarithmic_pdf (n, 0.4);
1845 test_lognormal (void)
1847 return gsl_ran_lognormal (r_global, 2.7, 1.3);
1851 test_lognormal_pdf (double x)
1853 return gsl_ran_lognormal_pdf (x, 2.7, 1.3);
1857 test_multinomial (void)
1860 const unsigned int sum_n = BINS;
1862 /* Test use of weights instead of probabilities. */
1863 const double p[] = { 2., 7., 1.};
1865 gsl_ran_multinomial ( r_global, K, sum_n, p, n);
1871 test_multinomial_pdf (unsigned int n_0)
1873 /* The margional distribution of just 1 variate is binomial. */
1875 /* Test use of weights instead of probabilities */
1876 double p[] = { 0.4, 1.6};
1877 const unsigned int sum_n = BINS;
1883 return gsl_ran_multinomial_pdf (K, p, n);
1888 test_multinomial_large (void)
1890 const unsigned int sum_n = BINS;
1891 unsigned int n[MULTI_DIM];
1892 const double p[MULTI_DIM] = { 0.2, 0.20, 0.17, 0.14, 0.12,
1893 0.07, 0.05, 0.04, 0.01, 0.00 };
1895 gsl_ran_multinomial ( r_global, MULTI_DIM, sum_n, p, n);
1901 test_multinomial_large_pdf (unsigned int n_0)
1903 return test_multinomial_pdf(n_0);
1907 test_negative_binomial (void)
1909 return gsl_ran_negative_binomial (r_global, 0.3, 20.0);
1913 test_negative_binomial_pdf (unsigned int n)
1915 return gsl_ran_negative_binomial_pdf (n, 0.3, 20.0);
1921 return gsl_ran_pascal (r_global, 0.8, 3);
1925 test_pascal_pdf (unsigned int n)
1927 return gsl_ran_pascal_pdf (n, 0.8, 3);
1934 return gsl_ran_pareto (r_global, 1.9, 2.75);
1938 test_pareto_pdf (double x)
1940 return gsl_ran_pareto_pdf (x, 1.9, 2.75);
1944 test_rayleigh (void)
1946 return gsl_ran_rayleigh (r_global, 1.9);
1950 test_rayleigh_pdf (double x)
1952 return gsl_ran_rayleigh_pdf (x, 1.9);
1956 test_rayleigh_tail (void)
1958 return gsl_ran_rayleigh_tail (r_global, 2.7, 1.9);
1962 test_rayleigh_tail_pdf (double x)
1964 return gsl_ran_rayleigh_tail_pdf (x, 2.7, 1.9);
1971 return gsl_ran_poisson (r_global, 5.0);
1975 test_poisson_pdf (unsigned int n)
1977 return gsl_ran_poisson_pdf (n, 5.0);
1981 test_poisson_large (void)
1983 return gsl_ran_poisson (r_global, 30.0);
1987 test_poisson_large_pdf (unsigned int n)
1989 return gsl_ran_poisson_pdf (n, 30.0);
1996 return gsl_ran_tdist (r_global, 1.75);
2000 test_tdist1_pdf (double x)
2002 return gsl_ran_tdist_pdf (x, 1.75);
2008 return gsl_ran_tdist (r_global, 12.75);
2012 test_tdist2_pdf (double x)
2014 return gsl_ran_tdist_pdf (x, 12.75);
2021 return gsl_ran_laplace (r_global, 2.75);
2025 test_laplace_pdf (double x)
2027 return gsl_ran_laplace_pdf (x, 2.75);
2033 return gsl_ran_weibull (r_global, 3.14, 2.75);
2037 test_weibull_pdf (double x)
2039 return gsl_ran_weibull_pdf (x, 3.14, 2.75);
2044 test_weibull1 (void)
2046 return gsl_ran_weibull (r_global, 2.97, 1.0);
2050 test_weibull1_pdf (double x)
2052 return gsl_ran_weibull_pdf (x, 2.97, 1.0);