5 #include <gnumeric-config.h>
7 #include <gnm-random.h>
11 #include <glib/gstdio.h>
20 * mathfunc.c: Mathematical functions.
23 * Morten Welinder <terra@gnome.org>
24 * Makoto Matsumoto and Takuji Nishimura (Mersenne Twister, see note in code)
25 * James Theiler. (See note 1.)
26 * Brian Gough. (See note 1.)
29 * NOTE 1: most of the random distribution code comes from the GNU Scientific
30 * Library (GSL), notably version 1.1.1. GSL is distributed under GPL licence,
31 * see COPYING. The relevant parts are copyright (C) 1996, 1997, 1998, 1999,
32 * 2000 James Theiler and Brian Gough.
37 /* ------------------------------------------------------------------------- */
38 /* http://www.math.keio.ac.jp/matumoto/CODES/MT2002/mt19937ar.c */
39 /* Imported by hand -- MW. */
42 A C-program for MT19937, with initialization improved 2002/1/26.
43 Coded by Takuji Nishimura and Makoto Matsumoto.
45 Before using, initialize the state by using init_genrand(seed)
46 or init_by_array(init_key, key_length).
48 Copyright (C) 1997 - 2002, Makoto Matsumoto and Takuji Nishimura,
51 Redistribution and use in source and binary forms, with or without
52 modification, are permitted provided that the following conditions
55 1. Redistributions of source code must retain the above copyright
56 notice, this list of conditions and the following disclaimer.
58 2. Redistributions in binary form must reproduce the above copyright
59 notice, this list of conditions and the following disclaimer in the
60 documentation and/or other materials provided with the distribution.
62 3. The names of its contributors may not be used to endorse or promote
63 products derived from this software without specific prior written
66 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
67 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
68 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
69 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
70 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
71 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
72 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
73 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
74 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
75 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
76 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
79 Any feedback is very welcome.
80 http://www.math.keio.ac.jp/matumoto/emt.html
81 email: matumoto@math.keio.ac.jp
88 /* Period parameters */
91 #define MATRIX_A 0x9908b0dfUL /* constant vector a */
92 #define UPPER_MASK 0x80000000UL /* most significant w-r bits */
93 #define LOWER_MASK 0x7fffffffUL /* least significant r bits */
95 static unsigned long mt
[N
]; /* the array for the state vector */
96 static int mti
=N
+1; /* mti==N+1 means mt[N] is not initialized */
98 /* initializes mt[N] with a seed */
99 static void init_genrand(unsigned long s
)
101 mt
[0]= s
& 0xffffffffUL
;
102 for (mti
=1; mti
<N
; mti
++) {
104 (1812433253UL * (mt
[mti
-1] ^ (mt
[mti
-1] >> 30)) + mti
);
105 /* See Knuth TAOCP Vol2. 3rd Ed. P.106 for multiplier. */
106 /* In the previous versions, MSBs of the seed affect */
107 /* only MSBs of the array mt[]. */
108 /* 2002/01/09 modified by Makoto Matsumoto */
109 mt
[mti
] &= 0xffffffffUL
;
110 /* for >32 bit machines */
114 /* initialize by an array with array-length */
115 /* init_key is the array for initializing keys */
116 /* key_length is its length */
117 static void mt_init_by_array(unsigned long init_key
[], int key_length
)
120 init_genrand(19650218UL);
122 k
= (N
>key_length
? N
: key_length
);
124 mt
[i
] = (mt
[i
] ^ ((mt
[i
-1] ^ (mt
[i
-1] >> 30)) * 1664525UL))
125 + init_key
[j
] + j
; /* non linear */
126 mt
[i
] &= 0xffffffffUL
; /* for WORDSIZE > 32 machines */
128 if (i
>=N
) { mt
[0] = mt
[N
-1]; i
=1; }
129 if (j
>=key_length
) j
=0;
131 for (k
=N
-1; k
; k
--) {
132 mt
[i
] = (mt
[i
] ^ ((mt
[i
-1] ^ (mt
[i
-1] >> 30)) * 1566083941UL))
133 - i
; /* non linear */
134 mt
[i
] &= 0xffffffffUL
; /* for WORDSIZE > 32 machines */
136 if (i
>=N
) { mt
[0] = mt
[N
-1]; i
=1; }
139 mt
[0] = 0x80000000UL
; /* MSB is 1; assuring non-zero initial array */
142 /* generates a random number on [0,0xffffffff]-interval */
143 static unsigned long genrand_int32(void)
146 static unsigned long mag01
[2]={0x0UL
, MATRIX_A
};
147 /* mag01[x] = x * MATRIX_A for x=0,1 */
149 if (mti
>= N
) { /* generate N words at one time */
152 if (mti
== N
+1) /* if init_genrand() has not been called, */
153 init_genrand(5489UL); /* a default initial seed is used */
155 for (kk
=0;kk
<N
-M
;kk
++) {
156 y
= (mt
[kk
]&UPPER_MASK
)|(mt
[kk
+1]&LOWER_MASK
);
157 mt
[kk
] = mt
[kk
+M
] ^ (y
>> 1) ^ mag01
[y
& 0x1UL
];
160 y
= (mt
[kk
]&UPPER_MASK
)|(mt
[kk
+1]&LOWER_MASK
);
161 mt
[kk
] = mt
[kk
+(M
-N
)] ^ (y
>> 1) ^ mag01
[y
& 0x1UL
];
163 y
= (mt
[N
-1]&UPPER_MASK
)|(mt
[0]&LOWER_MASK
);
164 mt
[N
-1] = mt
[M
-1] ^ (y
>> 1) ^ mag01
[y
& 0x1UL
];
173 y
^= (y
<< 7) & 0x9d2c5680UL
;
174 y
^= (y
<< 15) & 0xefc60000UL
;
181 /* generates a random number on [0,0x7fffffff]-interval */
182 long genrand_int31(void)
184 return (long)(genrand_int32()>>1);
187 /* generates a random number on [0,1]-real-interval */
188 double genrand_real1(void)
190 return genrand_int32()*(1.0/4294967295.0);
191 /* divided by 2^32-1 */
194 /* generates a random number on [0,1)-real-interval */
195 double genrand_real2(void)
197 return genrand_int32()*(1.0/4294967296.0);
198 /* divided by 2^32 */
201 /* generates a random number on (0,1)-real-interval */
202 double genrand_real3(void)
204 return (((double)genrand_int32()) + 0.5)*(1.0/4294967296.0);
205 /* divided by 2^32 */
208 /* generates a random number on [0,1) with 53-bit resolution*/
209 static double genrand_res53(void)
211 unsigned long a
=genrand_int32()>>5, b
=genrand_int32()>>6;
212 return(a
*67108864.0+b
)*(1.0/9007199254740992.0);
214 /* These real versions are due to Isaku Wada, 2002/01/09 added */
221 unsigned long init
[4]={0x123, 0x234, 0x345, 0x456}, length
=4;
222 init_by_array(init
, length
);
223 g_printerr("1000 outputs of genrand_int32()\n");
224 for (i
=0; i
<1000; i
++) {
225 g_printerr("%10lu ", genrand_int32());
226 if (i
%5==4) g_printerr("\n");
228 g_printerr("\n1000 outputs of genrand_real2()\n");
229 for (i
=0; i
<1000; i
++) {
230 g_printerr("%10.8f ", genrand_real2());
231 if (i
%5==4) g_printerr("\n");
244 /* ------------------------------------------------------------------------ */
247 mt_setup_seed (const char *seed
)
249 int len
= strlen (seed
);
251 unsigned long *longs
= g_new (unsigned long, len
+ 1);
253 /* We drop only one character into each long. */
254 for (i
= 0; i
< len
; i
++)
255 longs
[i
] = (unsigned char)seed
[i
];
256 mt_init_by_array (longs
, len
);
263 mt_setup_win32 (void)
265 /* See http://msdn.microsoft.com/en-us/library/Aa387694 */
266 typedef BOOLEAN (CALLBACK
* LPFNRTLGENRANDOM
) (void*,ULONG
);
267 LPFNRTLGENRANDOM MyRtlGenRandom
;
268 unsigned long buffer
[256];
270 gboolean res
= FALSE
;
272 hmod
= GetModuleHandle ("ADVAPI32.DLL");
276 MyRtlGenRandom
= (LPFNRTLGENRANDOM
)
277 GetProcAddress(hmod
, "SystemFunction036");
278 if (MyRtlGenRandom
&&
279 MyRtlGenRandom (buffer
, sizeof (buffer
))) {
280 mt_init_by_array (buffer
, G_N_ELEMENTS (buffer
));
291 /* ------------------------------------------------------------------------ */
294 random_01_mersenne (void)
296 size_t N
= (sizeof (gnm_float
) + sizeof (guint32
) - 1) / sizeof (guint32
);
303 for (n
= 0; n
< N
; n
++)
304 res
= (res
+ genrand_int32()) / 4294967296.0;
306 * It is conceivable that rounding turned the result
307 * into 1, so repeat in that case.
314 /* ------------------------------------------------------------------------ */
317 random_01_data (const unsigned char *data
)
322 for (ui
= 0; ui
< sizeof (gnm_float
); ui
++)
323 res
= (res
+ data
[ui
]) / 256;
327 /* ------------------------------------------------------------------------ */
329 static FILE *random_device_file
= NULL
;
330 #define RANDOM_DEVICE "/dev/urandom"
333 random_01_device (void)
335 static size_t bytes_left
= 0;
336 static unsigned char data
[32 * sizeof (gnm_float
)];
338 while (bytes_left
< sizeof (gnm_float
)) {
339 gssize items
= fread (data
+ bytes_left
, 1,
340 sizeof (data
) - bytes_left
,
343 g_warning ("Reading from %s failed; reverting to pseudo-random.",
345 return random_01_mersenne ();
350 bytes_left
-= sizeof (gnm_float
);
351 return random_01_data (data
+ bytes_left
);
354 /* ------------------------------------------------------------------------ */
362 static RandomSource random_src
= RS_UNDETERMINED
;
365 random_01_determine (void)
367 char const *seed
= g_getenv ("GNUMERIC_PRNG_SEED");
369 mt_setup_seed (seed
);
370 g_warning ("Using pseudo-random numbers.");
371 random_src
= RS_MERSENNE
;
376 if (mt_setup_win32 ()) {
377 random_src
= RS_MERSENNE
;
382 random_device_file
= g_fopen (RANDOM_DEVICE
, "rb");
383 if (random_device_file
) {
384 random_src
= RS_DEVICE
;
389 g_warning ("Using pseudo-random numbers.");
390 random_src
= RS_MERSENNE
;
397 if (random_src
== RS_UNDETERMINED
)
398 random_01_determine ();
400 switch (random_src
) {
401 case RS_UNDETERMINED
:
403 g_assert_not_reached ();
405 return random_01_mersenne ();
407 return random_01_device ();
411 /* ------------------------------------------------------------------------ */
416 * Returns: a N(0,1) distributed random number.
421 static gboolean has_saved
= FALSE
;
422 static gnm_float saved
;
428 gnm_float u
, v
, r2
, rsq
;
430 u
= 2 * random_01 () - 1;
431 v
= 2 * random_01 () - 1;
433 } while (r2
> 1 || r2
== 0);
435 rsq
= gnm_sqrt (-2 * gnm_log (r2
) / r2
);
445 random_lognormal (gnm_float zeta
, gnm_float sigma
)
447 return gnm_exp (sigma
* random_normal () + zeta
);
451 random_gaussian (gnm_float sigma
)
453 return sigma
* random_normal ();
457 * Generate a poisson distributed number.
460 random_poisson (gnm_float lambda
)
463 * This may not be optimal code, but it sure is easy to
464 * understand compared to R's code.
466 return qpois (random_01 (), lambda
, TRUE
, FALSE
);
470 * Generate a binomial distributed number.
473 random_binomial (gnm_float p
, gnm_float trials
)
475 return qbinom (random_01 (), trials
, p
, TRUE
, FALSE
);
479 * Generate a negative binomial distributed number.
482 random_negbinom (gnm_float p
, gnm_float f
)
484 return qnbinom (random_01 (), f
, p
, TRUE
, FALSE
);
488 * Generate an exponential distributed number.
491 random_exponential (gnm_float b
)
493 return -b
* gnm_log (random_01 ());
497 * Generate a bernoulli distributed number.
500 random_bernoulli (gnm_float p
)
502 gnm_float r
= random_01 ();
504 return (r
<= p
) ? 1.0 : 0.0;
508 * Generate a cauchy distributed number. From the GNU Scientific library 1.1.1.
509 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
512 random_cauchy (gnm_float a
)
520 return a
* gnm_tan (M_PIgnum
* u
);
524 * Generate a Weibull distributed number. From the GNU Scientific
526 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
529 random_weibull (gnm_float a
, gnm_float b
)
537 z
= gnm_pow (-gnm_log (x
), 1 / b
);
543 * Generate a Laplace (two-sided exponential probability) distributed number.
544 * From the GNU Scientific library 1.1.1.
545 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
548 random_laplace (gnm_float a
)
553 u
= 2 * random_01 () - 1.0;
557 return a
* gnm_log (-u
);
559 return -a
* gnm_log (u
);
563 random_laplace_pdf (gnm_float x
, gnm_float a
)
565 return (1 / (2 * a
)) * gnm_exp (-gnm_abs (x
) / a
);
569 * Generate a Rayleigh distributed number. From the GNU Scientific library
571 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
574 random_rayleigh (gnm_float sigma
)
582 return sigma
* gnm_sqrt (-2.0 * gnm_log (u
));
586 * Generate a Rayleigh tail distributed number. From the GNU Scientific library
587 * 1.1.1. The Rayleigh tail distribution has the form
588 * p(x) dx = (x / sigma^2) exp((a^2 - x^2)/(2 sigma^2)) dx
590 * for x = a ... +infty
592 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
595 random_rayleigh_tail (gnm_float a
, gnm_float sigma
)
603 return gnm_sqrt (a
* a
- 2.0 * sigma
* sigma
* gnm_log (u
));
606 /* The Gamma distribution of order a>0 is defined by:
608 * p(x) dx = {1 / \Gamma(a) b^a } x^{a-1} e^{-x/b} dx
610 * for x>0. If X and Y are independent gamma-distributed random
611 * variables of order a1 and a2 with the same scale parameter b, then
612 * X+Y has gamma distribution of order a1+a2.
614 * The algorithms below are from Knuth, vol 2, 2nd ed, p. 129.
618 gamma_frac (gnm_float a
)
620 /* This is exercise 16 from Knuth; see page 135, and the solution is
624 gnm_float p
= M_Egnum
/ (a
+ M_Egnum
);
627 gnm_float u
= random_01 ();
633 x
= gnm_pow (v
, 1 / a
);
637 q
= gnm_pow (x
, a
- 1);
639 } while (random_01 () >= q
);
645 gamma_large (gnm_float a
)
648 * Works only if a > 1, and is most efficient if a is large
650 * This algorithm, reported in Knuth, is attributed to Ahrens. A
651 * faster one, we are told, can be found in: J. H. Ahrens and
652 * U. Dieter, Computing 12 (1974) 223-246.
655 gnm_float sqa
, x
, y
, v
;
656 sqa
= gnm_sqrt (2 * a
- 1);
659 y
= gnm_tan (M_PIgnum
* random_01 ());
663 } while (v
> (1 + y
* y
) * gnm_exp ((a
- 1) * gnm_log (x
/ (a
- 1)) -
670 ran_gamma_int (gnm_float a
)
678 ua
= (unsigned int)a
;
680 for (i
= 0; i
< ua
; i
++)
681 prod
*= random_01 ();
684 * This handles the 0-probability event of getting
685 * an actual zero as well as the possibility of
690 return -gnm_log (prod
);
692 return gamma_large (a
);
696 * Generate a Gamma distributed number. From the GNU Scientific library 1.1.1.
697 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
700 random_gamma (gnm_float a
, gnm_float b
)
704 if (gnm_isnan (a
) || gnm_isnan (b
) || a
<= 0 || b
<= 0)
710 return b
* ran_gamma_int (na
);
712 return b
* gamma_frac (a
);
714 return b
* (ran_gamma_int (na
) + gamma_frac (a
- na
));
718 * Generate a Pareto distributed number. From the GNU Scientific library 1.1.1.
719 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
722 random_pareto (gnm_float a
, gnm_float b
)
730 return b
* gnm_pow (x
, -1 / a
);
734 * Generate a F-distributed number. From the GNU Scientific library 1.1.1.
735 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
738 random_fdist (gnm_float nu1
, gnm_float nu2
)
740 gnm_float Y1
= random_gamma (nu1
/ 2, 2.0);
741 gnm_float Y2
= random_gamma (nu2
/ 2, 2.0);
743 return (Y1
* nu2
) / (Y2
* nu1
);
747 * Generate a Beta-distributed number. From the GNU Scientific library 1.1.1.
748 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
751 random_beta (gnm_float a
, gnm_float b
)
753 gnm_float x1
= random_gamma (a
, 1.0);
754 gnm_float x2
= random_gamma (b
, 1.0);
756 return x1
/ (x1
+ x2
);
760 * Generate a Chi-Square-distributed number. From the GNU Scientific library
762 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
765 random_chisq (gnm_float nu
)
767 return 2 * random_gamma (nu
/ 2, 1.0);
771 * Generate a logistic-distributed number. From the GNU Scientific library
773 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
776 random_logistic (gnm_float a
)
782 } while (x
== 0 || x
== 1);
784 return a
* gnm_log (x
/ (1 - x
));
788 * Generate a geometric-distributed number. From the GNU Scientific library
790 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
793 random_geometric (gnm_float p
)
804 * Change from gsl version: we have support {0,1,2,...}
806 return gnm_floor (gnm_log (u
) / gnm_log1p (-p
));
810 random_hypergeometric (gnm_float n1
, gnm_float n2
, gnm_float t
)
812 return qhyper (random_01 (), n1
, n2
, t
, TRUE
, FALSE
);
817 * Generate a logarithmic-distributed number. From the GNU Scientific library
819 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
822 random_logarithmic (gnm_float p
)
839 q
= -gnm_expm1 (c
* u
);
842 return gnm_floor (1 + gnm_log (v
) / gnm_log (q
));
851 * Generate a T-distributed number. From the GNU Scientific library 1.1.1.
852 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
855 random_tdist (gnm_float nu
)
858 gnm_float Y1
= random_normal ();
859 gnm_float Y2
= random_chisq (nu
);
861 gnm_float t
= Y1
/ gnm_sqrt (Y2
/ nu
);
865 gnm_float Y1
, Y2
, Z
, t
;
867 Y1
= random_normal ();
868 Y2
= random_exponential (1 / (nu
/ 2 - 1));
870 Z
= Y1
* Y1
/ (nu
- 2);
871 } while (1 - Z
< 0 || gnm_exp (-Y2
- Z
) > (1 - Z
));
873 /* Note that there is a typo in Knuth's formula, the line below
874 * is taken from the original paper of Marsaglia, Mathematics
875 * of Computation, 34 (1980), p 234-256. */
877 t
= Y1
/ gnm_sqrt ((1 - 2 / nu
) * (1 - Z
));
883 * Generate a Type I Gumbel-distributed random number. From the GNU
884 * Scientific library 1.1.1.
885 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
888 random_gumbel1 (gnm_float a
, gnm_float b
)
896 return (gnm_log (b
) - gnm_log (-gnm_log (x
))) / a
;
900 * Generate a Type II Gumbel-distributed random number. From the GNU
901 * Scientific library 1.1.1.
902 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
905 random_gumbel2 (gnm_float a
, gnm_float b
)
913 return gnm_pow (-b
/ gnm_log (x
), 1 / a
);
917 * Generate a stable Levy-distributed random number. From the GNU
918 * Scientific library 1.1.1.
919 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough.
921 * The stable Levy probability distributions have the form
923 * p(x) dx = (1/(2 pi)) \int dt exp(- it x - |c t|^alpha)
925 * with 0 < alpha <= 2.
927 * For alpha = 1, we get the Cauchy distribution
928 * For alpha = 2, we get the Gaussian distribution with sigma = sqrt(2) c.
930 * Fromn Chapter 5 of Bratley, Fox and Schrage "A Guide to
931 * Simulation". The original reference given there is,
933 * J.M. Chambers, C.L. Mallows and B. W. Stuck. "A method for
934 * simulating stable random variates". Journal of the American
935 * Statistical Association, JASA 71 340-344 (1976).
938 random_levy (gnm_float c
, gnm_float alpha
)
940 gnm_float u
, v
, t
, s
;
946 u
= M_PIgnum
* (u
- 0.5);
948 if (alpha
== 1) { /* cauchy case */
954 v
= random_exponential (1.0);
957 if (alpha
== 2) { /* gaussian case */
958 t
= 2 * gnm_sin (u
) * gnm_sqrt (v
);
964 t
= gnm_sin (alpha
* u
) / gnm_pow (gnm_cos (u
), 1 / alpha
);
965 s
= gnm_pow (gnm_cos ((1 - alpha
) * u
) / v
, (1 - alpha
) / alpha
);
971 * The following routine for the skew-symmetric case was provided by
974 * The stable Levy probability distributions have the form
978 * = int dt exp(mu*i*t-|sigma*t|^alpha*(1-i*beta*sign(t)*tan(pi*alpha/2))) for
980 * = int dt exp(mu*i*t-|sigma*t|^alpha*(1+i*beta*sign(t)*2/pi*log(|t|))) for
983 * with 0<alpha<=2, -1<=beta<=1, sigma>0.
985 * For beta=0, sigma=c, mu=0, we get gsl_ran_levy above.
987 * For alpha = 1, beta=0, we get the Lorentz distribution
988 * For alpha = 2, beta=0, we get the Gaussian distribution
990 * See A. Weron and R. Weron: Computer simulation of Lévy alpha-stable
991 * variables and processes, preprint Technical University of Wroclaw.
992 * http://www.im.pwr.wroc.pl/~hugo/Publications.html
995 random_levy_skew (gnm_float c
, gnm_float alpha
, gnm_float beta
)
999 if (beta
== 0) /* symmetric case */
1000 return random_levy (c
, alpha
);
1006 V
= M_PIgnum
* (V
- 0.5);
1009 W
= random_exponential (1.0);
1013 X
= ((M_PI_2gnum
+ beta
* V
) * gnm_tan (V
) -
1014 beta
* gnm_log (M_PI_2gnum
* W
* gnm_cos (V
) /
1015 (M_PI_2gnum
+ beta
* V
))) / M_PI_2gnum
;
1016 return c
* (X
+ beta
* gnm_log (c
) / M_PI_2gnum
);
1018 gnm_float t
= beta
* gnm_tan (M_PI_2gnum
* alpha
);
1019 gnm_float B
= gnm_atan (t
) / alpha
;
1020 gnm_float S
= pow1p (t
* t
, 1 / (2 * alpha
));
1022 X
= S
* gnm_sin (alpha
* (V
+ B
)) / gnm_pow (gnm_cos (V
),
1024 * gnm_pow (gnm_cos (V
- alpha
* (V
+ B
)) / W
,
1025 (1 - alpha
) / alpha
);
1031 random_exppow_pdf (gnm_float x
, gnm_float a
, gnm_float b
)
1033 gnm_float lngamma
= lgamma1p (1 / b
);
1035 return (1 / (2 * a
)) * gnm_exp (-gnm_pow (gnm_abs (x
/ a
), b
) - lngamma
);
1039 * The exponential power probability distribution is
1041 * p(x) dx = (1/(2 a Gamma(1+1/b))) * exp(-|x/a|^b) dx
1043 * for -infty < x < infty. For b = 1 it reduces to the Laplace
1046 * The exponential power distribution is related to the gamma
1047 * distribution by E = a * pow(G(1/b),1/b), where E is an exponential
1048 * power variate and G is a gamma variate.
1050 * We use this relation for b < 1. For b >=1 we use rejection methods
1051 * based on the laplace and gaussian distributions which should be
1054 * See P. R. Tadikamalla, "Random Sampling from the Exponential Power
1055 * Distribution", Journal of the American Statistical Association,
1056 * September 1980, Volume 75, Number 371, pages 683-686.
1058 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
1062 random_exppow (gnm_float a
, gnm_float b
)
1064 /* See http://www.mcgill.ca/files/economics/propertiesandestimation.pdf */
1065 if (!(a
> 0) || gnm_isnan (b
))
1069 gnm_float u
= random_01 ();
1070 gnm_float v
= random_gamma (1 / b
, 1.0);
1071 gnm_float z
= a
* gnm_pow (v
, 1 / b
) ;
1078 return random_laplace (a
); /* Laplace distribution */
1080 /* Use laplace distribution for rejection method */
1081 gnm_float x
, y
, h
, ratio
, u
;
1083 /* Scale factor chosen by upper bound on ratio at b = 2 */
1084 gnm_float s
= 1.4489;
1086 x
= random_laplace (a
);
1087 y
= random_laplace_pdf (x
, a
);
1088 h
= random_exppow_pdf (x
, a
, b
);
1089 ratio
= h
/ (s
* y
);
1091 } while (u
> ratio
);
1094 } else if (b
== 2) /* Gaussian distribution */
1095 return random_gaussian (a
/ gnm_sqrt (2.0));
1097 /* Use gaussian for rejection method */
1098 gnm_float x
, y
, h
, ratio
, u
;
1099 const gnm_float sigma
= a
/ gnm_sqrt (2.0);
1101 /* Scale factor chosen by upper bound on ratio at b = infinity.
1102 * This could be improved by using a rational function
1103 * approximation to the bounding curve. */
1105 gnm_float s
= 2.4091 ; /* this is sqrt(pi) e / 2 */
1108 x
= random_gaussian (sigma
) ;
1109 y
= dnorm (x
, 0.0, gnm_abs (sigma
), FALSE
) ;
1110 h
= random_exppow_pdf (x
, a
, b
) ;
1111 ratio
= h
/ (s
* y
) ;
1113 } while (u
> ratio
);
1120 * Generate a Gaussian tail-distributed random number. From the GNU
1121 * Scientific library 1.1.1.
1122 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
1125 random_gaussian_tail (gnm_float a
, gnm_float sigma
)
1128 * Returns a gaussian random variable larger than a
1129 * This implementation does one-sided upper-tailed deviates.
1132 gnm_float s
= a
/ sigma
;
1135 /* For small s, use a direct rejection method. The limit s < 1
1136 * can be adjusted to optimise the overall efficiency */
1141 x
= random_gaussian (1.0);
1145 /* Use the "supertail" deviates from the last two steps
1146 * of Marsaglia's rectangle-wedge-tail method, as described
1147 * in Knuth, v2, 3rd ed, pp 123-128. (See also exercise 11,
1148 * p139, and the solution, p586.)
1158 x
= gnm_sqrt (s
* s
- 2 * gnm_log (v
));
1159 } while (x
* u
> s
);
1165 * Generate a Landau-distributed random number. From the GNU Scientific
1168 * Copyright (C) 2001 David Morrison
1170 * Adapted from the CERN library routines DENLAN, RANLAN, and DISLAN
1171 * as described in http://consult.cern.ch/shortwrups/g110/top.html.
1172 * Original author: K.S. K\"olbig.
1174 * The distribution is given by the complex path integral,
1176 * p(x) = (1/(2 pi i)) \int_{c-i\inf}^{c+i\inf} ds exp(s log(s) + x s)
1178 * which can be converted into a real integral over [0,+\inf]
1180 * p(x) = (1/pi) \int_0^\inf dt \exp(-t log(t) - x t) sin(pi t)
1184 random_landau (void)
1186 static gnm_float F
[983] = {
1188 * Add empty element [0] to account for difference
1189 * between C and Fortran convention for lower bound.
1191 00.000000, 00.000000, 00.000000, 00.000000, 00.000000,
1192 -2.244733, -2.204365, -2.168163, -2.135219, -2.104898,
1193 -2.076740, -2.050397, -2.025605, -2.002150, -1.979866,
1194 -1.958612, -1.938275, -1.918760, -1.899984, -1.881879,
1195 -1.864385, -1.847451, -1.831030, -1.815083, -1.799574,
1196 -1.784473, -1.769751, -1.755383, -1.741346, -1.727620,
1197 -1.714187, -1.701029, -1.688130, -1.675477, -1.663057,
1198 -1.650858, -1.638868, -1.627078, -1.615477, -1.604058,
1199 -1.592811, -1.581729, -1.570806, -1.560034, -1.549407,
1200 -1.538919, -1.528565, -1.518339, -1.508237, -1.498254,
1201 -1.488386, -1.478628, -1.468976, -1.459428, -1.449979,
1202 -1.440626, -1.431365, -1.422195, -1.413111, -1.404112,
1203 -1.395194, -1.386356, -1.377594, -1.368906, -1.360291,
1204 -1.351746, -1.343269, -1.334859, -1.326512, -1.318229,
1205 -1.310006, -1.301843, -1.293737, -1.285688, -1.277693,
1206 -1.269752, -1.261863, -1.254024, -1.246235, -1.238494,
1207 -1.230800, -1.223153, -1.215550, -1.207990, -1.200474,
1208 -1.192999, -1.185566, -1.178172, -1.170817, -1.163500,
1209 -1.156220, -1.148977, -1.141770, -1.134598, -1.127459,
1210 -1.120354, -1.113282, -1.106242, -1.099233, -1.092255,
1211 -1.085306, -1.078388, -1.071498, -1.064636, -1.057802,
1212 -1.050996, -1.044215, -1.037461, -1.030733, -1.024029,
1213 -1.017350, -1.010695, -1.004064, -0.997456, -0.990871,
1214 -0.984308, -0.977767, -0.971247, -0.964749, -0.958271,
1215 -0.951813, -0.945375, -0.938957, -0.932558, -0.926178,
1216 -0.919816, -0.913472, -0.907146, -0.900838, -0.894547,
1217 -0.888272, -0.882014, -0.875773, -0.869547, -0.863337,
1218 -0.857142, -0.850963, -0.844798, -0.838648, -0.832512,
1219 -0.826390, -0.820282, -0.814187, -0.808106, -0.802038,
1220 -0.795982, -0.789940, -0.783909, -0.777891, -0.771884,
1221 -0.765889, -0.759906, -0.753934, -0.747973, -0.742023,
1222 -0.736084, -0.730155, -0.724237, -0.718328, -0.712429,
1223 -0.706541, -0.700661, -0.694791, -0.688931, -0.683079,
1224 -0.677236, -0.671402, -0.665576, -0.659759, -0.653950,
1225 -0.648149, -0.642356, -0.636570, -0.630793, -0.625022,
1226 -0.619259, -0.613503, -0.607754, -0.602012, -0.596276,
1227 -0.590548, -0.584825, -0.579109, -0.573399, -0.567695,
1228 -0.561997, -0.556305, -0.550618, -0.544937, -0.539262,
1229 -0.533592, -0.527926, -0.522266, -0.516611, -0.510961,
1230 -0.505315, -0.499674, -0.494037, -0.488405, -0.482777,
1231 -0.477153, -0.471533, -0.465917, -0.460305, -0.454697,
1232 -0.449092, -0.443491, -0.437893, -0.432299, -0.426707,
1233 -0.421119, -0.415534, -0.409951, -0.404372, -0.398795,
1234 -0.393221, -0.387649, -0.382080, -0.376513, -0.370949,
1235 -0.365387, -0.359826, -0.354268, -0.348712, -0.343157,
1236 -0.337604, -0.332053, -0.326503, -0.320955, -0.315408,
1237 -0.309863, -0.304318, -0.298775, -0.293233, -0.287692,
1238 -0.282152, -0.276613, -0.271074, -0.265536, -0.259999,
1239 -0.254462, -0.248926, -0.243389, -0.237854, -0.232318,
1240 -0.226783, -0.221247, -0.215712, -0.210176, -0.204641,
1241 -0.199105, -0.193568, -0.188032, -0.182495, -0.176957,
1242 -0.171419, -0.165880, -0.160341, -0.154800, -0.149259,
1243 -0.143717, -0.138173, -0.132629, -0.127083, -0.121537,
1244 -0.115989, -0.110439, -0.104889, -0.099336, -0.093782,
1245 -0.088227, -0.082670, -0.077111, -0.071550, -0.065987,
1246 -0.060423, -0.054856, -0.049288, -0.043717, -0.038144,
1247 -0.032569, -0.026991, -0.021411, -0.015828, -0.010243,
1248 -0.004656, 00.000934, 00.006527, 00.012123, 00.017722,
1249 00.023323, 00.028928, 00.034535, 00.040146, 00.045759,
1250 00.051376, 00.056997, 00.062620, 00.068247, 00.073877,
1251 00.079511, 00.085149, 00.090790, 00.096435, 00.102083,
1252 00.107736, 00.113392, 00.119052, 00.124716, 00.130385,
1253 00.136057, 00.141734, 00.147414, 00.153100, 00.158789,
1254 00.164483, 00.170181, 00.175884, 00.181592, 00.187304,
1255 00.193021, 00.198743, 00.204469, 00.210201, 00.215937,
1256 00.221678, 00.227425, 00.233177, 00.238933, 00.244696,
1257 00.250463, 00.256236, 00.262014, 00.267798, 00.273587,
1258 00.279382, 00.285183, 00.290989, 00.296801, 00.302619,
1259 00.308443, 00.314273, 00.320109, 00.325951, 00.331799,
1260 00.337654, 00.343515, 00.349382, 00.355255, 00.361135,
1261 00.367022, 00.372915, 00.378815, 00.384721, 00.390634,
1262 00.396554, 00.402481, 00.408415, 00.414356, 00.420304,
1263 00.426260, 00.432222, 00.438192, 00.444169, 00.450153,
1264 00.456145, 00.462144, 00.468151, 00.474166, 00.480188,
1265 00.486218, 00.492256, 00.498302, 00.504356, 00.510418,
1266 00.516488, 00.522566, 00.528653, 00.534747, 00.540850,
1267 00.546962, 00.553082, 00.559210, 00.565347, 00.571493,
1268 00.577648, 00.583811, 00.589983, 00.596164, 00.602355,
1269 00.608554, 00.614762, 00.620980, 00.627207, 00.633444,
1270 00.639689, 00.645945, 00.652210, 00.658484, 00.664768,
1271 00.671062, 00.677366, 00.683680, 00.690004, 00.696338,
1272 00.702682, 00.709036, 00.715400, 00.721775, 00.728160,
1273 00.734556, 00.740963, 00.747379, 00.753807, 00.760246,
1274 00.766695, 00.773155, 00.779627, 00.786109, 00.792603,
1275 00.799107, 00.805624, 00.812151, 00.818690, 00.825241,
1276 00.831803, 00.838377, 00.844962, 00.851560, 00.858170,
1277 00.864791, 00.871425, 00.878071, 00.884729, 00.891399,
1278 00.898082, 00.904778, 00.911486, 00.918206, 00.924940,
1279 00.931686, 00.938446, 00.945218, 00.952003, 00.958802,
1280 00.965614, 00.972439, 00.979278, 00.986130, 00.992996,
1281 00.999875, 01.006769, 01.013676, 01.020597, 01.027533,
1282 01.034482, 01.041446, 01.048424, 01.055417, 01.062424,
1283 01.069446, 01.076482, 01.083534, 01.090600, 01.097681,
1284 01.104778, 01.111889, 01.119016, 01.126159, 01.133316,
1285 01.140490, 01.147679, 01.154884, 01.162105, 01.169342,
1286 01.176595, 01.183864, 01.191149, 01.198451, 01.205770,
1287 01.213105, 01.220457, 01.227826, 01.235211, 01.242614,
1288 01.250034, 01.257471, 01.264926, 01.272398, 01.279888,
1289 01.287395, 01.294921, 01.302464, 01.310026, 01.317605,
1290 01.325203, 01.332819, 01.340454, 01.348108, 01.355780,
1291 01.363472, 01.371182, 01.378912, 01.386660, 01.394429,
1292 01.402216, 01.410024, 01.417851, 01.425698, 01.433565,
1293 01.441453, 01.449360, 01.457288, 01.465237, 01.473206,
1294 01.481196, 01.489208, 01.497240, 01.505293, 01.513368,
1295 01.521465, 01.529583, 01.537723, 01.545885, 01.554068,
1296 01.562275, 01.570503, 01.578754, 01.587028, 01.595325,
1297 01.603644, 01.611987, 01.620353, 01.628743, 01.637156,
1298 01.645593, 01.654053, 01.662538, 01.671047, 01.679581,
1299 01.688139, 01.696721, 01.705329, 01.713961, 01.722619,
1300 01.731303, 01.740011, 01.748746, 01.757506, 01.766293,
1301 01.775106, 01.783945, 01.792810, 01.801703, 01.810623,
1302 01.819569, 01.828543, 01.837545, 01.846574, 01.855631,
1303 01.864717, 01.873830, 01.882972, 01.892143, 01.901343,
1304 01.910572, 01.919830, 01.929117, 01.938434, 01.947781,
1305 01.957158, 01.966566, 01.976004, 01.985473, 01.994972,
1306 02.004503, 02.014065, 02.023659, 02.033285, 02.042943,
1307 02.052633, 02.062355, 02.072110, 02.081899, 02.091720,
1308 02.101575, 02.111464, 02.121386, 02.131343, 02.141334,
1309 02.151360, 02.161421, 02.171517, 02.181648, 02.191815,
1310 02.202018, 02.212257, 02.222533, 02.232845, 02.243195,
1311 02.253582, 02.264006, 02.274468, 02.284968, 02.295507,
1312 02.306084, 02.316701, 02.327356, 02.338051, 02.348786,
1313 02.359562, 02.370377, 02.381234, 02.392131, 02.403070,
1314 02.414051, 02.425073, 02.436138, 02.447246, 02.458397,
1315 02.469591, 02.480828, 02.492110, 02.503436, 02.514807,
1316 02.526222, 02.537684, 02.549190, 02.560743, 02.572343,
1317 02.583989, 02.595682, 02.607423, 02.619212, 02.631050,
1318 02.642936, 02.654871, 02.666855, 02.678890, 02.690975,
1319 02.703110, 02.715297, 02.727535, 02.739825, 02.752168,
1320 02.764563, 02.777012, 02.789514, 02.802070, 02.814681,
1321 02.827347, 02.840069, 02.852846, 02.865680, 02.878570,
1322 02.891518, 02.904524, 02.917588, 02.930712, 02.943894,
1323 02.957136, 02.970439, 02.983802, 02.997227, 03.010714,
1324 03.024263, 03.037875, 03.051551, 03.065290, 03.079095,
1325 03.092965, 03.106900, 03.120902, 03.134971, 03.149107,
1326 03.163312, 03.177585, 03.191928, 03.206340, 03.220824,
1327 03.235378, 03.250005, 03.264704, 03.279477, 03.294323,
1328 03.309244, 03.324240, 03.339312, 03.354461, 03.369687,
1329 03.384992, 03.400375, 03.415838, 03.431381, 03.447005,
1330 03.462711, 03.478500, 03.494372, 03.510328, 03.526370,
1331 03.542497, 03.558711, 03.575012, 03.591402, 03.607881,
1332 03.624450, 03.641111, 03.657863, 03.674708, 03.691646,
1333 03.708680, 03.725809, 03.743034, 03.760357, 03.777779,
1334 03.795300, 03.812921, 03.830645, 03.848470, 03.866400,
1335 03.884434, 03.902574, 03.920821, 03.939176, 03.957640,
1336 03.976215, 03.994901, 04.013699, 04.032612, 04.051639,
1337 04.070783, 04.090045, 04.109425, 04.128925, 04.148547,
1338 04.168292, 04.188160, 04.208154, 04.228275, 04.248524,
1339 04.268903, 04.289413, 04.310056, 04.330832, 04.351745,
1340 04.372794, 04.393982, 04.415310, 04.436781, 04.458395,
1341 04.480154, 04.502060, 04.524114, 04.546319, 04.568676,
1342 04.591187, 04.613854, 04.636678, 04.659662, 04.682807,
1343 04.706116, 04.729590, 04.753231, 04.777041, 04.801024,
1344 04.825179, 04.849511, 04.874020, 04.898710, 04.923582,
1345 04.948639, 04.973883, 04.999316, 05.024942, 05.050761,
1346 05.076778, 05.102993, 05.129411, 05.156034, 05.182864,
1347 05.209903, 05.237156, 05.264625, 05.292312, 05.320220,
1348 05.348354, 05.376714, 05.405306, 05.434131, 05.463193,
1349 05.492496, 05.522042, 05.551836, 05.581880, 05.612178,
1350 05.642734, 05.673552, 05.704634, 05.735986, 05.767610,
1351 05.799512, 05.831694, 05.864161, 05.896918, 05.929968,
1352 05.963316, 05.996967, 06.030925, 06.065194, 06.099780,
1353 06.134687, 06.169921, 06.205486, 06.241387, 06.277630,
1354 06.314220, 06.351163, 06.388465, 06.426130, 06.464166,
1355 06.502578, 06.541371, 06.580553, 06.620130, 06.660109,
1356 06.700495, 06.741297, 06.782520, 06.824173, 06.866262,
1357 06.908795, 06.951780, 06.995225, 07.039137, 07.083525,
1358 07.128398, 07.173764, 07.219632, 07.266011, 07.312910,
1359 07.360339, 07.408308, 07.456827, 07.505905, 07.555554,
1360 07.605785, 07.656608, 07.708035, 07.760077, 07.812747,
1361 07.866057, 07.920019, 07.974647, 08.029953, 08.085952,
1362 08.142657, 08.200083, 08.258245, 08.317158, 08.376837,
1363 08.437300, 08.498562, 08.560641, 08.623554, 08.687319,
1364 08.751955, 08.817481, 08.883916, 08.951282, 09.019600,
1365 09.088889, 09.159174, 09.230477, 09.302822, 09.376233,
1366 09.450735, 09.526355, 09.603118, 09.681054, 09.760191,
1367 09.840558, 09.922186, 10.005107, 10.089353, 10.174959,
1368 10.261958, 10.350389, 10.440287, 10.531693, 10.624646,
1369 10.719188, 10.815362, 10.913214, 11.012789, 11.114137,
1370 11.217307, 11.322352, 11.429325, 11.538283, 11.649285,
1371 11.762390, 11.877664, 11.995170, 12.114979, 12.237161,
1372 12.361791, 12.488946, 12.618708, 12.751161, 12.886394,
1373 13.024498, 13.165570, 13.309711, 13.457026, 13.607625,
1374 13.761625, 13.919145, 14.080314, 14.245263, 14.414134,
1375 14.587072, 14.764233, 14.945778, 15.131877, 15.322712,
1376 15.518470, 15.719353, 15.925570, 16.137345, 16.354912,
1377 16.578520, 16.808433, 17.044929, 17.288305, 17.538873,
1378 17.796967, 18.062943, 18.337176, 18.620068, 18.912049,
1379 19.213574, 19.525133, 19.847249, 20.180480, 20.525429,
1380 20.882738, 21.253102, 21.637266, 22.036036, 22.450278,
1381 22.880933, 23.329017, 23.795634, 24.281981, 24.789364,
1382 25.319207, 25.873062, 26.452634, 27.059789, 27.696581,
1383 28.365274, 29.068370, 29.808638, 30.589157, 31.413354,
1384 32.285060, 33.208568, 34.188705, 35.230920, 36.341388,
1385 37.527131, 38.796172, 40.157721, 41.622399, 43.202525,
1386 44.912465, 46.769077, 48.792279, 51.005773, 53.437996,
1387 56.123356, 59.103894
1389 gnm_float X
, U
, V
, RANLAN
;
1399 if (I
>= 70 && I
<= 800)
1400 RANLAN
= F
[I
] + U
* (F
[I
+ 1] - F
[I
]);
1401 else if (I
>= 7 && I
<= 980)
1402 RANLAN
= F
[I
] + U
* (F
[I
+ 1] - F
[I
] -
1403 0.25 * (1 - U
) * (F
[I
+ 2] - F
[I
+ 1] -
1408 RANLAN
= ((0.99858950 + (3.45213058E1
+ 1.70854528E1
* U
) * U
) /
1409 (1 + (3.41760202E1
+ 4.01244582 * U
) * U
)) *
1410 ( -gnm_log ( -0.91893853 - V
) - 1);
1415 RANLAN
= (1.00060006 + 2.63991156E2
*
1416 U
+ 4.37320068E3
* V
) /
1417 ((1 + 2.57368075E2
* U
+ 3.41448018E3
* V
) * U
);
1419 RANLAN
= (1.00001538 + 6.07514119E3
* U
+
1421 ((1 + 6.06511919E3
* U
+ 6.94021044E5
* V
) * U
);
1428 /* ------------------------------------------------------------------------ */
1431 * Generate a skew-normal distributed random number.
1433 * based on the information provided at
1434 * http://azzalini.stat.unipd.it/SN/faq-r.html
1439 random_skew_normal (gnm_float a
)
1442 gnm_float delta
= a
/ gnm_sqrt(1 + a
* a
);
1443 gnm_float u
= random_normal ();
1444 gnm_float v
= random_normal ();
1446 result
= delta
* u
+ gnm_sqrt (1-delta
*delta
) * v
;
1448 return ((u
< 0.) ? -result
: result
);
1452 /* ------------------------------------------------------------------------ */
1455 * Generate a skew-t distributed random number.
1457 * based on the information provided at
1458 * http://azzalini.stat.unipd.it/SN/faq-r.html
1463 random_skew_tdist (gnm_float nu
, gnm_float a
)
1465 gnm_float chi
= random_chisq (nu
);
1466 gnm_float z
= random_skew_normal (a
);
1468 return (z
* gnm_sqrt(nu
/chi
));
1471 /* ------------------------------------------------------------------------ */