random.c

贝叶斯算法：盲分离技术

    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 + 0.5 * v * v;

        p = 0.0;   // arbitrary
        r = (ycap - y1) + m1 * x1;
        if( r <= 0.0 )               // Rslope(0, ycap, -d/2, -u-a*v, v*v/2)
            q = m0 - sqrt(2.0 * d * height);
        else                         // Rslope(0, ycap, -1/2., -u, 0)
            q = (v - u) + 2.0 * height / (v - sqrt(2.0 * ycap));
    }
    else                             // apex in (0,x1)
    {
        xcap = m0 / d;
        ycap = height + 0.5 * (d * xcap * xcap + v * v);

        r = (ycap - y0) - m0 * xcap;
        if( r <= 0.0 )               // Lslope(xcap, ycap, -d/2, -u-a*v, v*v/2)
            p = sqrt(2.0 * d * height);
        else
            p = (ycap - y0) / xcap;
        r = (ycap - y1) - m1 * (xcap - x1);
        if( r <= 0.0 )               // Rslope(xcap, ycap, -d/2, -u-a*v, v*v/2)
            q = -sqrt(2.0 * d * height);
        else                         // Rslope(xcap, ycap, -1/2., -u, 0)
            q = (a * (u - v) + (1.0 + a) * v) * a / d
                 - sqrt(2.0 * r + (xcap - x1) * (xcap - x1));
    }

    *xc = xcap;
    *p0 = p;
    *q0 = q;
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  Positive2D
// 
// Purpose:   Set bi-exponential sampler for Positive2 in case a < 0, v >= 0
//      x = inf       |
//                    | exp  - x^2 / 2 - u x
//      x = 0         |
//                 
// History:   JS         31 Dec 2001, 3 Oct 2002
//-----------------------------------------------------------------------------
// 
static void Positive2D(
double  u,        // I    >= v >= 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)
{

// Apex necessarily at 0 
    *xc = 0.0;
    *p0 = 0.0;   // arbitrary
    *q0 = -u - sqrt(2.0 * height);           // Rslope(0, height, -1/2., -u, 0)
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  Ran1posneg
// 
// Purpose:   Draw sample from normalised form of
//
//                   exp( - u|x| - f x - A x^2 / 2 )
//
// Method:
//   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.
//
// History:   JS         31 Dec 2001, 3 Oct 2000
//-----------------------------------------------------------------------------
// 
double Ran1posneg( //   O  sample value
Rand_t  Rand,      // I O  random generator state
double  f,         // I    coeff of x
double  u,         // I    coeff of |x|,  > 0
double  A)         // I    coeff of x^2,  >= 0
{
static const double height = 0.5000;
    double  x, xcap, ycap;
    double  accept, p, q, r, m9, m0;

// Special case of no curvature
    if( A <= 0.0 )
    {
        if( Randouble(Rand) < 0.5 * (1.0 + f / u) )
            return log(Randouble(Rand)) / (u - f);
        else
            return -log(Randouble(Rand)) / (u + f);
    }
// Scale general case to unit curvature
    f /= sqrt(A);
    u /= sqrt(A);
    m9 = u - f;
    m0 = -u - f;
// Set dominating bi-exponential sampler
    if( m9 < 0.0 )          // Apex in x -ve
    {
        xcap = m9;
        ycap = height + 0.5 * xcap * xcap;
        p = sqrt(2.0 * height);             // Lslope(xcap, ycap, -1/2., m9, 0)
        if( ycap <= m9 * xcap )
            q = -sqrt(2.0 * height);        // Rslope(xcap, ycap, -1/2., m9, 0)
        else
        {
            r = ycap - m0 * xcap;
            if( r <= 0.0 )
                q = ycap / xcap;
            else                            // Rslope(xcap, ycap, -1/2., m0, 0)
                q = -2.0 * u - sqrt(2.0 * r + m9 * m9);
        }
    }
    else if( m0 <= 0.0 )   // Apex at x = 0
    {
        xcap = 0.0;
        ycap = height;
        p = m9 + sqrt(2.0 * height);        // Lslope(xcap, ycap, -1/2., m9, 0)
        q = m0 - sqrt(2.0 * height);        // Rslope(xcap, ycap, -1/2., m0, 0)
    }
    else                   // Apex in x +ve
    {
        xcap = m0;
        ycap = height + 0.5 * xcap * xcap;
        if( ycap <= m0 * xcap )
            p = sqrt(2.0 * height);         // Lslope(xcap, ycap, -1/2., m0, 0)
        else
        {
            r = ycap - m9 * xcap;
            if( r <= 0.0 )
                p = ycap / xcap;
            else                            // Lslope(xcap, ycap, -1/2., m9, 0)
                p = 2.0 * u + sqrt(2.0 * r + m0 * m0);
        }
        q = -sqrt(2.0 * height);            // Rslope(xcap, ycap, -1/2., m0, 0)
    }
// Sample x
    do
    {
// accept = -log(sampler)
        if( Randouble(Rand) < p / (p - q) )
        {
            x = xcap + log(Randouble(Rand)) / q;
            accept = -(ycap + q * (x - xcap));
        }
        else
        {
            x = xcap + log(Randouble(Rand)) / p;
            accept = -(ycap + p * (x - xcap));
        }
// accept += log(true)  
        if( x > 0.0 )
            accept += x * (- u - f - 0.5 * x);
        else
            accept += x * (u - f - 0.5 * x);
    } while( log(Randouble(Rand)) > accept );
// Exit
    return x / sqrt(A);
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  Ran2posneg
// 
// Purpose:   Draw sample from normalised form of
//
//    exp( - u|x| - v|y| - f x - g y - (A11 x^2 + 2 A12 x y + A22 y^2) / 2 )
//
// History:   JS         31 Dec 2001, 3 Oct 2002
//-----------------------------------------------------------------------------
// 
void Ran2posneg(
Rand_t  Rand,     // I O  random generator state
double  f,        // I    coeff of x
double  g,        // I    coeff of y
double  u,        // I    coeff of |x|,  > 0
double  v,        // 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 = Ran1posneg(Rand, f, u, A11);
        *y = Ran1posneg(Rand, g, v, A22);
        return;
    }
    f /= sqrt(A11);
    u /= sqrt(A11);
    g /= sqrt(A22);
    v /= sqrt(A22);
    a = A12 / sqrt(A11 * A22);
    if( a > 1.0 )
        a = 1.0;
    if( a < -1.0 )
        a = -1.0;

    if( f <= 0.0 )
    {
        if( g <= 0.0 )
        {
            if( -g - v >= -f - u )
                Posneg2(Rand, -f, -g, u, v, -a, x, y);
            else
                Posneg2(Rand, -g, -f, v, u, -a, y, x);
        }
        else
        {
            if( g - v >= -f - u )
                Posneg2(Rand, -f, g, u, v, a, x, y);
            else
                Posneg2(Rand, g, -f, v, u, a, y, x);
            *y = - *y;
        }
    }
    else
    {
        if( g <= 0.0 )
        {
            if( -g - v >= f - u )
                Posneg2(Rand, f, -g, u, v, a, x, y);
            else
                Posneg2(Rand, -g, f, v, u, a, y, x);
        }
        else
        {
            if( g - v >= f - u )
                Posneg2(Rand, f, g, u, v, -a, x, y);
            else
                Posneg2(Rand, g, f, v, u, -a, y, x);
            *y = - *y;
        }
        *x = - *x;
    }

    *x /= sqrt(A11);
    *y /= sqrt(A22);
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  Posneg2
// 
// Purpose:   Draw sample from normalised form of
//
//        exp( - u|x| - v|y| + f x + g y - (x^2 - 2 a x y + y^2) / 2 )
//
//            with  f >= 0, g >= 0, f - u <= g - v
//            as well as u > 0, v > 0, 0 < |a| <= 1
//
// Method:
//   Apart from the simple special case  u - f > 0.3 or so, v - g > 0.3 or so:
//
//                 | exp -x^2/2 + f x - u|x| + (g - v + a x)^2/2 ,  g-v+ax > 0
//   Sample x from | exp -x^2/2 + f x - u|x|                     ,  |g+ax| < v
//                 | exp -x^2/2 + f x - u|x| + (g + v + a x)^2/2 ,  g+v+ax < 0
//   then
//                 | exp -(y - g + v - a x)^2 / 2                ,  g-v+ax > 0
//   Sample y from | exp - y^2 / 2                               ,  |g+ax| < v
//                 | exp -(y - g - v - a x)^2 / 2                ,  g+v+ax < 0
//   with
//                 | exp - v (|y| - y)                           ,  g-v+ax > 0
//   Acceptance    | exp - v |y| + (g + a x) y                   ,  |g+ax| < v
//                 | exp - v (|y| + y)                           ,  g+v+ax < 0
//
//   Provided  f-u <= g-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 Posneg2(
Rand_t  Rand,     // I O  random generator state
double  f,        // I    coeff of x,    >= 0
double  g,        // I    coeff of y,    >= 0
double  u,        // I    coeff of |x|,  > 0
double  v,        // I    coeff of |y|,  > 0
double  a,        // I    coeff of xy,   -1 < 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 - g > 0.3000 )
// Special case: sample around (0,0) from exp  f x + g y - u|x| - v|y|
    {
        do
        {
            x = (Randouble(Rand) > 0.5 * (1.0 - f / u))
               ? log(Randouble(Rand)) / (f - u)
               : log(Randouble(Rand)) / (f + u);
            y = (Randouble(Rand) > 0.5 * (1.0 - g / v))
               ? log(Randouble(Rand)) / (g - v)
               : log(Randouble(Rand)) / (g + v);
            accept = a * x * y - 0.5 * (x * x + y * y);
        } while( accept < log(Randouble(Rand)) );
        *x0 = x;
        *y0 = y;
        return;
    }

// General bi-exponential for x
    do
    {
        if( a > 0.0 )
        {
            if( g > v )
                Posneg2A(f+u, f-u, g+v, g-v, a, height, &xcap, &p, &q);
            else
                Posneg2B(f+u, f-u, g+v, g-v, a, height, &xcap, &p, &q);
        }
        else
        {
            if( g > v )
                Posneg2C(f+u, f-u, g+v, g-v, a, height, &xcap, &p, &q);
            else
                Posneg2D(f+u, f-u, g+v, g-v, a, height, &xcap, &p, &q);
        }

// Sample x
        do
        {
            s = (Randouble(Rand) < p / (p - q)) ? q : p;
            dx = log(Randouble(Rand)) / s;
            x = dx + xcap;
            if( x > 0.0 )
                accept = dx * (f - u - 0.5 * (xcap + x));
            else
                accept = dx * (f + u - 0.5 * (xcap + x));
            r = g - v + a * x;
            if( r > 0.0 )
                accept += a * dx * (r - 0.5 * a * dx);
            r = g + v + a * x;
            if( r < 0.0 )
                accept += a * dx * (r - 0.5 * a * dx);
            accept -= height + s * dx;
        } while( log(Randouble(Rand)) > accept );

// Sample y|x
        if( g - v + a * x > 0.0 )
        {
            y = g - v + a * x + Rangauss(Rand);
            accept = (y >= 0.0) ? 0.0 : 2.0 * v * y;