random.c

贝叶斯算法：盲分离技术

                  & 0xfffff000) ^ 0x00000800;   // switch lowest (2^-53) bit ON
    return  ((double)hi + (double)lo * SHIFT32) * SHIFT32;
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  Rangauss
// 
// Purpose:   Sample from Gaussian N(0,1)
// 
// Notes: (1) Cannot overflow;  |value| < 10 standard deviations
//        (2) 2nd sample y*sqrt(...) could be available, best with y = ...+ 0.5
// 
// History:   JS          3 Oct 2002  Half of Box-Muller
//                       17 Dec 2002  Use Rand[3] for extra 40% speedup
//-----------------------------------------------------------------------------
// 
double Rangauss(       //   O  Value
Rand_t   Rand)         // I O  Random generator state
{
    static const double RR = 4611686018427387904.0;   // 2^62
    double x, y, r;
    do
    {
        x = Ranint(Rand) + 0.5;              // (-2^31,2^31), not 0
        y = (int)Rand[3];                    // [-2^31,2^31)
        r = x * x + y * y;
    } while( r >= RR );                      // disc, radius 2^31
    return  x * sqrt(2.0 * log(RR / r) / r); // Waste 2nd Box-Muller sample
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  Rancauchy
// 
// Purpose:   Draw sample from
//                prob(x)  =  1 / (1 + x * x)
// 
// Notes:     Cannot overflow;  |value| <= 2^53/pi
// 
// History:   JS         3 Oct 2002
//-----------------------------------------------------------------------------
// 
double  Rancauchy(     //   O  Value
Rand_t   Rand)         // I O  Random generator state
{
    double  t = tan(HalfPi * (Randouble(Rand) - 0.5));
    return (1.0 - t * t) / (2.0 * t);
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  Rangamma
// 
// Purpose:   Draw sample from
//                             -1+c  -x
//                prob(x)  =  x     e  / Gamma(c)   for x in [0, inf)
// 
// Notes:     x=0 can be returned, especially if c is small, but x=inf cannot.
// 
// History:   JS         24 Jan 1994, 19 Oct 1995, 24 Aug 1996, 3 Oct 2002
//-----------------------------------------------------------------------------
// 
double  Rangamma(      //   O  Value
Rand_t   Rand,         // I O  Random generator state
double   c)            // I    Exponent
{
    double a, q, r, g, e;

    g = c - 1.0;
    if (g > 0.0)
    {
        e = sqrt(g + c) / g;
        do
        {
            do
            {
                r = Rancauchy(Rand);
                q = e * r + 1.0;
            } while (q <= 0.0);
        } while( log(Randouble(Rand) / (1.0 + r * r))
                > g * (log(q) - q + 1.0) );             // Accept?
        q *= g;
    }
    else
    {
        r = 1.0 / (1.0 + c * exp(-1.0));                // prob(bounding p < 1)
        do
        {
            if( Randouble(Rand) >= r )                  // Bounding x > 1
            {
                q = 1.0 - log(Randouble(Rand));         // exp(-q) in (1,inf)
                a = g * log(q);                         // log(Accept ratio)
            }
            else                                        // Bounding x <= 1
            {
                q = exp(log(Randouble(Rand)) / c);      // q^(-1+c) in (0,1)
                a = -q;                                 // log(Accept ratio)
            }
        } while( log(Randouble(Rand)) > a );
    }
    return q;
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  Ranpoiss
// 
// Purpose:   Draw sample from Poisson
//
//                           -c  j
//                prob(j) = e   c / j!               ( 0 <= c <~ 4000000000 )
// 
// Notes:     Distribution is truncated at 2^32-1,
//            so input parameter c should not approach 2^32.
//
// History:   JS         15 May 1998, 24 Mar 2001, 3 Oct 2002
//-----------------------------------------------------------------------------
// 
unsigned   Ranpoiss(   //   O  Value
Rand_t   Rand,         // I O  Random generator state
double   c)            // I    Mean
{
    unsigned j, k;
    double   a, b, r, t, y;

    if( c <= 0.0 )
        return 0;
    if( c < 70.0 )                             // Direct is faster for small c
    {
        j = k = (unsigned)c;
        t = exp(j * log(c) - c - logfactorial(j));
        for( ; ; )
        {
            r = Randouble(Rand);
            if( t >= r )
                return j;
            a = b = t;
            do
            {
                if( a > b )
                {
                    if( (t += a *= j-- / c) >= r )
                        return j;
                }
                else
                {
                    if( (t += b *= c / ++k) >= r )
                        return k;
                }
            } while( b > 0.0 );            // Trap failure from rounding error
        }
    }
    else
    {
        b = sqrt(2.0 * c);
        do
        {
            do
            {
                y = Rancauchy(Rand);
                r = c + b * y;
            } while( r < 0.0 );
            j = (unsigned)r;
            t = j * log(c) - logfactorial(j) - c;
        } while( log(Randouble(Rand) / (1.731 * (1.0 + y * y) * b)) > t );
        return j;
    }
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  Ranbinom
// 
// Purpose:   Draw sample from binomial
//
//                             n!     j      n-j
//                prob(j) = -------- p  (1-p)              0 <= j <= n
//                          j!(n-j)!
// 
// History:   JS         24 Mar 2001, 3 Oct 2002
//-----------------------------------------------------------------------------
// 
unsigned  Ranbinom(    //   O  Value
Rand_t   Rand,         // I O  Random generator state
unsigned n,            // I    Range
double   p)            // I    Mean/Range
{
    unsigned j;
    double   a, b, g, q, r, t, u, w, y;

    if( n == 0 )
        return 0;
    if( p <= 0.0 )
        return 0;
    if( p >= 1.0 )
        return n;
    q = 1.0 - p;

    if( n < 100 )
    {
        u = w = j = (unsigned)(p * n + p);
        t = exp(j * log(p) + (n-j) * log(q)
                 + logfactorial(n) - logfactorial(j) - logfactorial(n-j));
        for( ; ; )
        {
            r = Randouble(Rand);
            if( t >= r )
                return j;
            a = b = t;
            do
            {
                if( a > b )
                {
                    a *= u / p;
                    u -= 1.0;
                    a *= q / (n - u);
                    t += a;
                    if( t >= r )
                        return (unsigned)u;
                }
                else
                {
                    b *= (n - w) / q;
                    w += 1.0;
                    b *= p / w;
                    t += b;
                    if( t >= r )
                        return (unsigned)w;
                }
            } while( a + b > 0.0 );      // Trap failure from rounding error
        }
    }
    else
    {
        b = sqrt(2.0 * n * p * q + 0.25);
        g = logfactorial(n);
        do
        {
            do
            {
                y = Rancauchy(Rand);
                r = n * p + b * y + 0.5;
            } while( r < 0.0 || r >= n + 1.0 );
            j = (unsigned)r;
            t = g - logfactorial(j) - logfactorial(n - j)
               + j * log(p) + (n - j) * log(q);
        } while( log(Randouble(Rand) / (0.8190 * (1.0 + y * y) * b)) > t );
        return j;
    }
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  Ranbeta
//
// Purpose:   Draw sample from beta distribution
//                             r       n-r
//                prob(x)  =  x (1 - x)   / Z ,     Z = (n+1)! / r! (n-r)!
//
// Method:    Rejection beneath Cauchy of phenomenologically adequate width.
// 
// History:   John Skilling   28 Nov 2000
//-----------------------------------------------------------------------------
// 
double Ranbeta(       //   O  Value
Rand_t   Rand,        // I O  Random generator state
int      r,           // I    Number of successes >= 0
int      n)           // I    Total number        >= r
{
    double p, q, x, y, c;
    x = (double)n;
    if( r <= 0 )
        x = 1.0 - exp(log(Randouble(Rand)) / (x + 1.0));
    else if( r >= n )
        x = exp(log(Randouble(Rand)) / (x + 1.0));
    else
    {    
        p = (double)r / x;
        q = 1.0 - p;
        c = x / (3.0 * p * q);             // less than curvature at top
        do
        {
            do
            {
                y = Rancauchy(Rand) / sqrt(c);
                x = p + y;                 // x from Cauchy upper bound 
            }
            while( x <= 0.0 || x >= 1.0 ); // x within range
            y = r * log(x / p) + (n-r) * log((1.0 - x) / q)   // log(true)
               + log(1.0 + c * y * y);                        // -log(Cauchy)
        } while( log(Randouble(Rand)) > y );     // accept?
    }
    return x;
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  Rangeom
// 
// Purpose:   Draw integer from truncated geometric distribution 
//                                              j
//                    prob(j)   proportional   a    for j in [0,N-1]
//
//                    where  a = alpha/(alpha+1),   alpha > 0.0
//
//            Special case N = 0 means N = 2^32 for which
//
//                    <j> = alpha ,   var(j) = alpha(alpha+1)
// 
// History:   John Skilling   3 Apr 2001  Linearly truncated version
//                           17 Dec 2002  Hard truncation
//-----------------------------------------------------------------------------
// 
unsigned Rangeom(     //   O  Value
Rand_t   Rand,        // I O  Random generator state
double   alpha,       // I    Parameter
unsigned N)           // I    Supremum (0 interpreted as 2^32)
{
    double   a, r;

    if( N == 1 || alpha <= 0.0 )
        return 0;
    if( alpha * DBL_EPSILON >= 0.5 )
        return Rangrid(Rand, N);

    a = alpha / (alpha + 1.0);
    do
    {
        r = Randouble(Rand);
        if( N )
            r *= 1.0 - pow(a, (double)N);
        r = log(1.0 - r) / log(a);
    } while( r >= N );            // safety
    return (int)r;
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  Rangrid
//
// Purpose:   Random integer from [0,bound-1].
//
// Method:    Development of Ranint.
//
// History:   John Skilling   6 May 1995 - 15 Jan 2002
//-----------------------------------------------------------------------------
// 
unsigned   Rangrid(    //   O  Value
Rand_t   Rand,         // I O  Random generator state
unsigned bound)        // I    Supremum (0 interpreted as 2^32)
{
    unsigned   i, M, N;
    if( bound == 0 )
        return (unsigned)Ranint(Rand);
    N = (unsigned)(-1) / bound;
    M = N * bound;
    do  i = (unsigned)Ranint(Rand);
    while( i >= M );
    return i / N;
 }

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  Ranperm
//
// Purpose:   Random permutation of {0,1,2,...,n-1}.
//
// History:   John Skilling   23 May 1998, 24 Jan 1999, 9 Nov 2000
//-----------------------------------------------------------------------------
// 
void Ranperm( 
Rand_t   Rand,        // I O  Random generator state
int      n,           // I    Dimension
int*     perm)        //   O  Output permutation
{
    int    i, m;
    for( m = 0; m < n; ++m )