📄 random.c
字号:
{
i = Rangrid(Rand, (unsigned)(m + 1));
if( i < m )
perm[m] = perm[i];
perm[i] = m;
}
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: Ran2gauss
//
// Purpose: Draw sample from normalised form of
//
// exp( - g.x - x.A.x / 2 ) ( A non-negative )
//
// History: JS 31 Dec 2001
//-----------------------------------------------------------------------------
//
void Ran2gauss(
Rand_t Rand, // I O Random generator state
double g1, // I x gradient at origin
double g2, // I y gradient at origin
double A11, // I xx curvature >= 0
double A12, // I xy curvature, A12*A12 <= A11*A22
double A22, // I yy curvature >= 0
double* x, // O x sample position
double* y) // O y sample position
{
double d = A11 * A22 - A12 * A12;
double xbar = (g2 * A12 - g1 * A22) / d;
double ybar = (g1 * A12 - g2 * A11) / d;
*x = Rangauss(Rand) * sqrt(A22 / d);
*y = Rangauss(Rand) / sqrt(A22) - *x * A12 / A22;
*x += xbar;
*y += ybar;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: Ran1pos
//
// Purpose: Draw positive sample (x +ve) from normalised form of
//
// exp( - g x - A x^2 / 2 )
//
// History: JS 30 Nov 1997
//-----------------------------------------------------------------------------
//
double Ran1pos( // O sample value
Rand_t Rand, // I O random generator state
double g, // I coeff of x
double A) // I coeff of x^2, >= 0
{
double a, x;
a = sqrt(A);
if( 0.3 * a > g ) // Reject from completed normal curve
{
do x = Rangauss(Rand) / a - g / A;
while( x < 0.0 );
}
else // Reject from bounding exponential
{
do x = -log(Randouble(Rand)) / g;
while( -log(Randouble(Rand)) < A * x * x / 2.0 );
}
return x;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: Ran2pos
//
// Purpose: Draw sample in positive quadrant from normalised form of
//
// exp( - g1 x - g2 y - (A11 x^2 + 2 A12 x y + A22 y^2) / 2 )
//
// History: JS 31 Dec 2001, 3 Oct 2002
//-----------------------------------------------------------------------------
//
void Ran2pos(
Rand_t Rand, // I O random generator state
double g1, // I coeff of |x|, > 0
double g2, // I coeff of |y|, > 0
double A11, // I coeff of x^2, >= 0
double A12, // I coeff of xy, |A12| < sqrt(A11*A22)
double A22, // I coeff of y^2, >= 0
double* x, // O sample
double* y) // O sample
{
double a;
if( A12 * A12 <= DBL_EPSILON * A11 * A22 )
{
*x = Ran1pos(Rand, g1, A11);
*y = Ran1pos(Rand, g2, A22);
return;
}
g1 /= sqrt(A11);
g2 /= sqrt(A22);
a = A12 / sqrt(A11 * A22);
if( a > 1.0 )
a = 1.0;
if( a < -1.0 )
a = -1.0;
if( g2 <= g1 )
Positive2(Rand, g1, g2, -a, x, y);
else
Positive2(Rand, g2, g1, -a, y, x);
*x /= sqrt(A11);
*y /= sqrt(A22);
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: Positive2
//
// Purpose: Draw sample in positive quadrant from normalised form of
//
// exp( - u x - v y - (x^2 - 2 a x y + y^2) / 2 )
//
// with u >= v as well as |a| < 1
//
// Method:
// Apart from the simple special case u > 0.3 or so, v > 0.3 or so:
//
// Sample x from | exp -x^2/2 - u x + (a x - v)^2/2 , ax > v
// | exp -x^2/2 - u x , ax < v
//
// then
// Sample y from | exp -(y - a x + v)^2 / 2 , ax > v
// | exp - y^2 / 2 , ax < v
// with
// Acceptance | 1 , ax > v
// | exp (a x - v)y , ax < v
//
// Provided u >= v, these algorithms always sample with O(1) efficiency.
//
// General sampler for ANY logconcave function of x (such as the above) is
// exp( height + ycap + q * (x - xcap) ), x > xcap (q < 0)
// exp( height + ycap + p * (x - xcap) ), x < xcap (p > 0)
// where (xcap,ycap) is apex of required curve over all x, height = O(1) +ve,
// and p,q define left and right exponentials tangent to the required curve
// A x^2 + B x + C (with A < 0), according to
// p = Lslope(x,y, A,B,C) = 2*A*x + B + 2 * sqrt(-A *(y - A*x*x - B*x - C))
// q = Rslope(x,y, A,B,C) = 2*A*x + B - 2 * sqrt(-A *(y - A*x*x - B*x - C))
// This method always samples with O(1) efficiency.
//
// Subsidiary codes crafted to evaluate all sqrt arguments definitely +ve,
// left slope p definitely +ve, and right slope q definitely -ve.
//
// History: JS 31 Dec 2001, 3 Oct 2002
//-----------------------------------------------------------------------------
//
static void Positive2(
Rand_t Rand, // I O random generator state
double u, // I coeff of x
double v, // I coeff of y
double a, // I coeff of xy, 0 < |a| <= 1
double* x0, // O sample
double* y0) // O sample
{
double height = 0.5000;
double x, y, accept, xcap, p, q, r, s, dx;
if( v > 0.3000 )
// Special case: sample around (0,0) from exp - u x - v y
{
do
{
x = -log(Randouble(Rand)) / u;
y = -log(Randouble(Rand)) / v;
accept = a * x * y - 0.5 * (x * x + y * y);
} while( accept < log(Randouble(Rand)) );
}
else
// General sampler ....
{
// Bi-exponential for x
for( ; ; )
{
if( a > 0.0 )
{
if( v < 0.0 )
Positive2A(u, v, a, height, &xcap, &p, &q);
else
Positive2B(u, v, a, height, &xcap, &p, &q);
}
else
{
if( v < 0.0 )
Positive2C(u, v, a, height, &xcap, &p, &q);
else
Positive2D(u, height, &xcap, &p, &q);
}
// Sample x
for( ; ; )
{
s = (xcap == 0.0 || Randouble(Rand) < p / (p - q)) ? q : p;
dx = log(Randouble(Rand)) / s;
x = dx + xcap;
if( x <= 0.0 )
continue;
accept = -dx * (u + xcap + 0.5 * dx);
r = a * x - v;
if( r > 0.0 )
accept += a * dx * (0.5 * a * (x + xcap) - v);
accept -= height + s * dx;
if( log(Randouble(Rand)) < accept )
break;
}
// Sample y|x
if( - v + a * x > 0.0 )
{
y = - v + a * x + Rangauss(Rand);
if( y >= 0.0 )
break;
}
else
{
y = Rangauss(Rand);
if( y >= 0.0 && log(Randouble(Rand)) < (- v + a * x) * y )
break;
}
}
}
// Exit
*x0 = x;
*y0 = y;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: Positive2A
//
// Purpose: Set bi-exponential sampler for Positive2 in case a > 0, v <= 0
// x = inf |
// | exp - x^2 / 2 - u x + (a x - v)^2 / 2
// x = 0 |
//
// History: JS 31 Dec 2001, 3 Oct 2002
//-----------------------------------------------------------------------------
//
static void Positive2A(
double u, // I >= v
double v, // I <= 0
double a, // I > 0
double height, // I height of sampler cusp above target maximum
double* xc, // O target maximum
double* p0, // O sampler exp( yc + p0 * (x-xc) ), x < xc (p0 > 0)
double* q0) // O sampler exp( yc + q0 * (x-xc) ), x > xc (q0 < 0)
{
double y0, m0, xcap, ycap;
double d, p, q;
d = 1.0 - a * a;
if( d < DBL_EPSILON )
d = DBL_EPSILON;
// Set boundaries (0,y0)...
// and slopes m0
y0 = 0.5 * v * v;
m0 = -(u + a * v);
// Find apex from slope=0
if( m0 <= 0.0 ) // apex at 0
{
xcap = 0.0;
ycap = height + y0;
p = 0.0; // arbitrary
q = m0 - sqrt(2.0 * d * height); // Rslope(xcap, ycap, -d/2, m0, v*v/2)
}
else // apex in (0,inf)
{
xcap = m0 / d;
ycap = height + 0.5 * (m0 * xcap + v * v);
if( ycap - y0 <= m0 * xcap )
p = sqrt(2.0 * d * height); // Lslope(xcap, ycap, -d/2, m0, v*v/2)
else
p = (ycap - y0) / xcap;
q = -sqrt(2.0 * d * height); // Rslope(xcap, ycap, -d/2, m0, v*v/2)
}
*xc = xcap;
*p0 = p;
*q0 = q;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: Positive2B
//
// Purpose: Set bi-exponential sampler for Positive2 in case a > 0, v >= 0
//
// | exp - x^2 / 2 - u x + (a x - v)^2 / 2
// x2 = v/a |
// | exp - x^2 / 2 - u x
// x = 0 |
//
// History: JS 31 Dec 2001, 3 Oct 2002
//-----------------------------------------------------------------------------
//
static void Positive2B(
double u, // I >= v
double v, // I >= 0
double a, // I > 0
double height, // I height of sampler cusp above target maximum
double* xc, // O target maximum
double* p0, // O sampler exp( yc + p0 * (x-xc) ), x < xc (p0 > 0)
double* q0) // O sampler exp( yc + q0 * (x-xc) ), x > xc (q0 < 0)
{
double x2, y2, m2;
double d, q, r;
d = 1.0 - a * a;
if( d < DBL_EPSILON )
d = DBL_EPSILON;
// Set boundaries (0,y0)...(x2,y2)...
// and slopes +ve,-ve m2
x2 = v / a;
y2 = x2 * (- u - 0.5 * x2);
m2 = - u - x2;
// Apex necessarily at 0
r = (height - y2) + m2 * x2;
if( r <= 0.0 ) // Rslope(0, height, -1/2., -u, 0)
q = -u - sqrt(2.0 * height);
else // Rslope(0, height, -d/2, -u-a*v, v*v/2)
q = -(u + a * v) - sqrt(2.0 * d * (d * height + a * a * r));
*xc = 0.0;
*p0 = 0.0; // arbitrary
*q0 = q;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: Positive2C
//
// Purpose: Set bi-exponential sampler for Positive2 in case a < 0, v <= 0
//
// | exp - x^2 / 2 - u x
// x1 = -v/-a |
// | exp - x^2 / 2 - u x + (a x - v)^2 / 2
// x = 0 |
//
// History: JS 31 Dec 2001, 3 Oct 2002
//-----------------------------------------------------------------------------
//
static void Positive2C(
double u, // I >= v
double v, // I <= 0
double a, // I < 0
double height, // I height of sampler cusp above target maximum
double* xc, // O target maximum
double* p0, // O sampler exp( yc + p0 * (x-xc) ), x < xc (p0 > 0)
double* q0) // O sampler exp( yc + q0 * (x-xc) ), x > xc (q0 < 0)
{
double x1, y0, y1, m0, m1, xcap, ycap;
double d, p, q, r;
d = 1.0 - a * a;
if( d < DBL_EPSILON )
d = DBL_EPSILON;
// Set boundaries (0,y0)...(x1,y1)...
// and slopes m0 m1
x1 = v / a;
y1 = x1 * (- u - 0.5 * x1);
m1 = - u - x1; // -ve
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -