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

   1 /* randist/levy.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 <math.h>
  22 #include <gsl/gsl_math.h>
  23 #include <gsl/gsl_rng.h>
  24 #include <gsl/gsl_randist.h>
  25
  26 /* The stable Levy probability distributions have the form
  27
  28    p(x) dx = (1/(2 pi)) \int dt exp(- it x - |c t|^alpha)
  29
  30    with 0 < alpha <= 2.
  31
  32    For alpha = 1, we get the Cauchy distribution
  33    For alpha = 2, we get the Gaussian distribution with sigma = sqrt(2) c.
  34
  35    Fromn Chapter 5 of Bratley, Fox and Schrage "A Guide to
  36    Simulation". The original reference given there is,
  37
  38    J.M. Chambers, C.L. Mallows and B. W. Stuck. "A method for
  39    simulating stable random variates". Journal of the American
  40    Statistical Association, JASA 71 340-344 (1976).
  41
  42    */
  43
  44 double
  45 gsl_ran_levy (const gsl_rng * r, const double c, const double alpha)
  46 {
  47   double u, v, t, s;
  48
  49   u = M_PI * (gsl_rng_uniform_pos (r) - 0.5);
  50
  51   if (alpha == 1)               /* cauchy case */
  52     {
  53       t = tan (u);
  54       return c * t;
  55     }
  56
  57   do
  58     {
  59       v = gsl_ran_exponential (r, 1.0);
  60     }
  61   while (v == 0);
  62
  63   if (alpha == 2)             /* gaussian case */
  64     {
  65       t = 2 * sin (u) * sqrt(v);
  66       return c * t;
  67     }
  68
  69   /* general case */
  70
  71   t = sin (alpha * u) / pow (cos (u), 1 / alpha);
  72   s = pow (cos ((1 - alpha) * u) / v, (1 - alpha) / alpha);
  73
  74   return c * t * s;
  75 }
  76
  77
  78 /* The following routine for the skew-symmetric case was provided by
  79    Keith Briggs.
  80
  81    The stable Levy probability distributions have the form
  82
  83    2*pi* p(x) dx
  84
  85      = \int dt exp(mu*i*t-|sigma*t|^alpha*(1-i*beta*sign(t)*tan(pi*alpha/2))) for alpha!=1
  86      = \int dt exp(mu*i*t-|sigma*t|^alpha*(1+i*beta*sign(t)*2/pi*log(|t|)))   for alpha==1
  87
  88    with 0<alpha<=2, -1<=beta<=1, sigma>0.
  89
  90    For beta=0, sigma=c, mu=0, we get gsl_ran_levy above.
  91
  92    For alpha = 1, beta=0, we get the Lorentz distribution
  93    For alpha = 2, beta=0, we get the Gaussian distribution
  94
  95    See A. Weron and R. Weron: Computer simulation of Lévy alpha-stable
  96    variables and processes, preprint Technical University of Wroclaw.
  97    http://www.im.pwr.wroc.pl/~hugo/Publications.html
  98
  99 */
 100
 101 double
 102 gsl_ran_levy_skew (const gsl_rng * r, const double c,
 103                    const double alpha, const double beta)
 104 {
 105   double V, W, X;
 106
 107   if (beta == 0)  /* symmetric case */
 108     {
 109       return gsl_ran_levy (r, c, alpha);
 110     }
 111
 112   V = M_PI * (gsl_rng_uniform_pos (r) - 0.5);
 113
 114   do
 115     {
 116       W = gsl_ran_exponential (r, 1.0);
 117     }
 118   while (W == 0);
 119
 120   if (alpha == 1)
 121     {
 122       X = ((M_PI_2 + beta * V) * tan (V) -
 123            beta * log (M_PI_2 * W * cos (V) / (M_PI_2 + beta * V))) / M_PI_2;
 124       return c * (X + beta * log (c) / M_PI_2);
 125     }
 126   else
 127     {
 128       double t = beta * tan (M_PI_2 * alpha);
 129       double B = atan (t) / alpha;
 130       double S = pow (1 + t * t, 1/(2 * alpha));
 131
 132       X = S * sin (alpha * (V + B)) / pow (cos (V), 1 / alpha)
 133         * pow (cos (V - alpha * (V + B)) / W, (1 - alpha) / alpha);
 134       return c * X;
 135     }
 136 }