|
Packit |
67cb25 |
/* randist/gauss.c
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Copyright (C) 1996, 1997, 1998, 1999, 2000, 2006, 2007 James Theiler, Brian Gough
|
|
Packit |
67cb25 |
* Copyright (C) 2006 Charles Karney
|
|
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 |
/* Of the two methods provided below, I think the Polar method is more
|
|
Packit |
67cb25 |
* efficient, but only when you are actually producing two random
|
|
Packit |
67cb25 |
* deviates. We don't produce two, because then we'd have to save one
|
|
Packit |
67cb25 |
* in a static variable for the next call, and that would screws up
|
|
Packit |
67cb25 |
* re-entrant or threaded code, so we only produce one. This makes
|
|
Packit |
67cb25 |
* the Ratio method suddenly more appealing.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* [Added by Charles Karney] We use Leva's implementation of the Ratio
|
|
Packit |
67cb25 |
* method which avoids calling log() nearly all the time and makes the
|
|
Packit |
67cb25 |
* Ratio method faster than the Polar method (when it produces just one
|
|
Packit |
67cb25 |
* result per call). Timing per call (gcc -O2 on 866MHz Pentium,
|
|
Packit |
67cb25 |
* average over 10^8 calls)
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* Polar method: 660 ns
|
|
Packit |
67cb25 |
* Ratio method: 368 ns
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
*/
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Polar (Box-Mueller) method; See Knuth v2, 3rd ed, p122 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_ran_gaussian (const gsl_rng * r, const double sigma)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double x, y, r2;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
do
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* choose x,y in uniform square (-1,-1) to (+1,+1) */
|
|
Packit |
67cb25 |
x = -1 + 2 * gsl_rng_uniform_pos (r);
|
|
Packit |
67cb25 |
y = -1 + 2 * gsl_rng_uniform_pos (r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* see if it is in the unit circle */
|
|
Packit |
67cb25 |
r2 = x * x + y * y;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
while (r2 > 1.0 || r2 == 0);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Box-Muller transform */
|
|
Packit |
67cb25 |
return sigma * y * sqrt (-2.0 * log (r2) / r2);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Ratio method (Kinderman-Monahan); see Knuth v2, 3rd ed, p130.
|
|
Packit |
67cb25 |
* K+M, ACM Trans Math Software 3 (1977) 257-260.
|
|
Packit |
67cb25 |
*
|
|
Packit |
67cb25 |
* [Added by Charles Karney] This is an implementation of Leva's
|
|
Packit |
67cb25 |
* modifications to the original K+M method; see:
|
|
Packit |
67cb25 |
* J. L. Leva, ACM Trans Math Software 18 (1992) 449-453 and 454-455. */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_ran_gaussian_ratio_method (const gsl_rng * r, const double sigma)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double u, v, x, y, Q;
|
|
Packit |
67cb25 |
const double s = 0.449871; /* Constants from Leva */
|
|
Packit |
67cb25 |
const double t = -0.386595;
|
|
Packit |
67cb25 |
const double a = 0.19600;
|
|
Packit |
67cb25 |
const double b = 0.25472;
|
|
Packit |
67cb25 |
const double r1 = 0.27597;
|
|
Packit |
67cb25 |
const double r2 = 0.27846;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
do /* This loop is executed 1.369 times on average */
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
/* Generate a point P = (u, v) uniform in a rectangle enclosing
|
|
Packit |
67cb25 |
the K+M region v^2 <= - 4 u^2 log(u). */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* u in (0, 1] to avoid singularity at u = 0 */
|
|
Packit |
67cb25 |
u = 1 - gsl_rng_uniform (r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* v is in the asymmetric interval [-0.5, 0.5). However v = -0.5
|
|
Packit |
67cb25 |
is rejected in the last part of the while clause. The
|
|
Packit |
67cb25 |
resulting normal deviate is strictly symmetric about 0
|
|
Packit |
67cb25 |
(provided that v is symmetric once v = -0.5 is excluded). */
|
|
Packit |
67cb25 |
v = gsl_rng_uniform (r) - 0.5;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Constant 1.7156 > sqrt(8/e) (for accuracy); but not by too
|
|
Packit |
67cb25 |
much (for efficiency). */
|
|
Packit |
67cb25 |
v *= 1.7156;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Compute Leva's quadratic form Q */
|
|
Packit |
67cb25 |
x = u - s;
|
|
Packit |
67cb25 |
y = fabs (v) - t;
|
|
Packit |
67cb25 |
Q = x * x + y * (a * y - b * x);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Accept P if Q < r1 (Leva) */
|
|
Packit |
67cb25 |
/* Reject P if Q > r2 (Leva) */
|
|
Packit |
67cb25 |
/* Accept if v^2 <= -4 u^2 log(u) (K+M) */
|
|
Packit |
67cb25 |
/* This final test is executed 0.012 times on average. */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
while (Q >= r1 && (Q > r2 || v * v > -4 * u * u * log (u)));
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return sigma * (v / u); /* Return slope */
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_ran_gaussian_pdf (const double x, const double sigma)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double u = x / fabs (sigma);
|
|
Packit |
67cb25 |
double p = (1 / (sqrt (2 * M_PI) * fabs (sigma))) * exp (-u * u / 2);
|
|
Packit |
67cb25 |
return p;
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_ran_ugaussian (const gsl_rng * r)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return gsl_ran_gaussian (r, 1.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_ran_ugaussian_ratio_method (const gsl_rng * r)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return gsl_ran_gaussian_ratio_method (r, 1.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
gsl_ran_ugaussian_pdf (const double x)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
return gsl_ran_gaussian_pdf (x, 1.0);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|