方法太多了,给你两个代码:
ratio method:
<br />
<br />
double Normal::gaussian_ratio_method (const double sigma)<br />
{<br />
double u, v, x, y, Q;<br />
const double s = 0.449871; /* Constants from Leva */<br />
const double t = -0.386595;<br />
const double a = 0.19600;<br />
const double b = 0.25472;<br />
const double r1 = 0.27597;<br />
const double r2 = 0.27846;<br />
<br />
do /* This loop is executed 1.369 times on average */<br />
{<br />
/* Generate a point P = (u, v) uniform in a rectangle enclosing<br />
the K+M region v^2 <= - 4 u^2 log(u). */<br />
<br />
/* u in (0, 1] to avoid singularity at u = 0 */<br />
u = 1 - rng_->next();<br />
<br />
/* v is in the asymmetric interval [-0.5, 0.5). However v = -0.5<br />
is rejected in the last part of the while clause. The<br />
resulting Normal deviate is strictly symmetric about 0<br />
(provided that v is symmetric once v = -0.5 is excluded). */<br />
v = rng_->next() - 0.5;<br />
<br />
/* Constant 1.7156 > sqrt(8/e) (for accuracy); but not by too<br />
much (for efficiency). */<br />
v *= 1.7156;<br />
<br />
/* Compute Leva's quadratic form Q */<br />
x = u - s;<br />
y = fabs (v) - t;<br />
Q = x * x + y * (a * y - b * x);<br />
<br />
/* Accept P if Q < r1 (Leva) */<br />
/* Reject P if Q > r2 (Leva) */<br />
/* Accept if v^2 <= -4 u^2 log(u) (K+M) */<br />
/* This final test is executed 0.012 times on average. */<br />
}<br />
while (Q >= r1 && (Q > r2 || v * v > -4 * u * u * log (u)));<br />
<br />
return sigma * (v / u); /* Return slope */<br />
}<br />
<br />
<br />
box-mueller method:
<br />
<br />
double Normal::gaussian_pola_method (const double sigma)<br />
{<br />
double x, y, r2;<br />
<br />
do<br />
{<br />
/* choose x,y in uniform square (-1,-1) to (+1,+1) */<br />
x = -1.0 + 2.0 * rng_->next_pos();<br />
y = -1.0 + 2.0 * rng_->next_pos();<br />
<br />
/* see if it is in the unit circle */<br />
r2 = x * x + y * y;<br />
}<br />
while (r2 > 1.0 || r2 == 0);<br />
<br />
/* Box-Muller transform */<br />
return sigma * y * sqrt (-2.0 * log (r2) / r2);<br />
}<br />
<br />
我用的是Ziggurat Method by George Marsaglia and Wai Wan Tsang
代码太长不贴了。