|
Packit |
67cb25 |
/* randist/levy.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is free software; you can redistribute it and/or modify
|
|
Packit |
67cb25 |
* it under the terms of the GNU General Public License as published by
|
|
Packit |
67cb25 |
* the Free Software Foundation; either version 3 of the License, or (at
|
|
Packit |
67cb25 |
* your option) any later version.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* This program is distributed in the hope that it will be useful, but
|
|
Packit |
67cb25 |
* WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
Packit |
67cb25 |
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
|
|
Packit |
67cb25 |
* General Public License for more details.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* You should have received a copy of the GNU General Public License
|
|
Packit |
67cb25 |
* along with this program; if not, write to the Free Software
|
|
Packit |
67cb25 |
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
#include <config.h>
|
|
Packit |
67cb25 |
#include <math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_rng.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_randist.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* The stable Levy probability distributions have the form
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
p(x) dx = (1/(2 pi)) \int dt exp(- it x - |c t|^alpha)
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with 0 < alpha <= 2.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
For alpha = 1, we get the Cauchy distribution
|
|
Packit |
67cb25 |
For alpha = 2, we get the Gaussian distribution with sigma = sqrt(2) c.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
Fromn Chapter 5 of Bratley, Fox and Schrage "A Guide to
|
|
Packit |
67cb25 |
Simulation". The original reference given there is,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
J.M. Chambers, C.L. Mallows and B. W. Stuck. "A method for
|
|
Packit |
67cb25 |
simulating stable random variates". Journal of the American
|
|
Packit |
67cb25 |
Statistical Association, JASA 71 340-344 (1976).
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_ran_levy (const gsl_rng * r, const double c, const double alpha)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double u, v, t, s;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
u = M_PI * (gsl_rng_uniform_pos (r) - 0.5);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (alpha == 1) /* cauchy case */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
t = tan (u);
|
|
Packit |
67cb25 |
return c * t;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
do
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
v = gsl_ran_exponential (r, 1.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
while (v == 0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (alpha == 2) /* gaussian case */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
t = 2 * sin (u) * sqrt(v);
|
|
Packit |
67cb25 |
return c * t;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* general case */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
t = sin (alpha * u) / pow (cos (u), 1 / alpha);
|
|
Packit |
67cb25 |
s = pow (cos ((1 - alpha) * u) / v, (1 - alpha) / alpha);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return c * t * s;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* The following routine for the skew-symmetric case was provided by
|
|
Packit |
67cb25 |
Keith Briggs.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
The stable Levy probability distributions have the form
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
2*pi* p(x) dx
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
= \int dt exp(mu*i*t-|sigma*t|^alpha*(1-i*beta*sign(t)*tan(pi*alpha/2))) for alpha!=1
|
|
Packit |
67cb25 |
= \int dt exp(mu*i*t-|sigma*t|^alpha*(1+i*beta*sign(t)*2/pi*log(|t|))) for alpha==1
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
with 0<alpha<=2, -1<=beta<=1, sigma>0.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
For beta=0, sigma=c, mu=0, we get gsl_ran_levy above.
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
For alpha = 1, beta=0, we get the Lorentz distribution
|
|
Packit |
67cb25 |
For alpha = 2, beta=0, we get the Gaussian distribution
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
See A. Weron and R. Weron: Computer simulation of Lévy alpha-stable
|
|
Packit |
67cb25 |
variables and processes, preprint Technical University of Wroclaw.
|
|
Packit |
67cb25 |
http://www.im.pwr.wroc.pl/~hugo/Publications.html
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_ran_levy_skew (const gsl_rng * r, const double c,
|
|
Packit |
67cb25 |
const double alpha, const double beta)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double V, W, X;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (beta == 0) /* symmetric case */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return gsl_ran_levy (r, c, alpha);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
V = M_PI * (gsl_rng_uniform_pos (r) - 0.5);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
do
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
W = gsl_ran_exponential (r, 1.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
while (W == 0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
if (alpha == 1)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
X = ((M_PI_2 + beta * V) * tan (V) -
|
|
Packit |
67cb25 |
beta * log (M_PI_2 * W * cos (V) / (M_PI_2 + beta * V))) / M_PI_2;
|
|
Packit |
67cb25 |
return c * (X + beta * log (c) / M_PI_2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
else
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double t = beta * tan (M_PI_2 * alpha);
|
|
Packit |
67cb25 |
double B = atan (t) / alpha;
|
|
Packit |
67cb25 |
double S = pow (1 + t * t, 1/(2 * alpha));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
X = S * sin (alpha * (V + B)) / pow (cos (V), 1 / alpha)
|
|
Packit |
67cb25 |
* pow (cos (V - alpha * (V + B)) / W, (1 - alpha) / alpha);
|
|
Packit |
67cb25 |
return c * X;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
}
|