bayesys3.c

贝叶斯算法：盲分离技术

    CALL( BayesInit(Common, Objects) )
// Initialise prior objects
    CALL( PriorInit(Common, Objects, Links, NEXTRA, Lextra) )

// Enter at prior (infinite temperature)
    dcool = BIG;
    ntable[0] = j = 0;
    ctable[0] = Ltable[0] = Evid = 0.0;
    do
    {
        Common->Nsystem ++;
// Update Ocean
        CALL( FillOcean(Common, Objects, Ocean, NOCEAN) )
// Controlled cooling
        cold = Common->cool;                       // previous value
        CALL( Control(Common, Objects, NEXTRA, Lextra, &t) )
        Common->cool += t;
        if( Common->cool > cold + 2.0000 * dcool ) // MAX allowable cooling to
            Common->cool = cold + 2.0000 * dcool;  // prevent fast acceleration
// User-limited cooling, at posterior (= 1) or beyond (> 1)
// Reset any nuisance parameters
// Collect statistics and display diagnostics
        CALL( UserMonitor(Common, Objects) )       // should not increase dcool
        finish = CALLvalue;                        // enables user to exit
        if( Common->cool > BIG )                   // avoid overflow
            Common->cool = BIG;
        dcool = Common->cool - cold;
// Cool ensemble by copying (user's nuisance parameters won't be copied)
        CALL( Anneal(Common, Objects, NEXTRA, dcool) )
// Evolve
        CALL( MCMCengines(Common, Objects, Links) )
// Evidence and Information diagnostics
        Lbar = 0.0;
        for( k = 0; k < ENSEMBLE; k++ )
            Lbar += Objects[k].Lhood;
        Lbar /= ENSEMBLE;
        if( j == NTABLE - 1 )
        {
            NTABLE += NTABLE / 2;
            REALLOC(Ltable, NTABLE, double)
            REALLOC(ctable, NTABLE, double)
            REALLOC(ntable, NTABLE, int)
        }
        if( Common->cool > ctable[j] )   // new table value
        {
            j++;
            ctable[j] = Common->cool;
            Ltable[j] = Lbar;
            ntable[j] = 1;
            Evid += (ctable[j] - ctable[j-1]) * Ltable[j];
        }
        else                             // update existing top value
        {
            Evid -= (ctable[j] - ctable[j-1]) * Ltable[j];
            Ltable[j] = (ntable[j] * Ltable[j] + Lbar) / (ntable[j] + 1);
            ntable[j] ++;
            Evid += (ctable[j] - ctable[j-1]) * Ltable[j];
        }
// Prune tables to satisfy known constraint: logL increasing function of cool
// If logL tries to decrease, feed the deficit backwards
        while( j > 1 && Ltable[j-1] >= Ltable[j] )
        {
            Evid -= (ctable[j] - ctable[j-1]) * Ltable[j];
            j--;
            Evid -= (ctable[j] - ctable[j-1]) * Ltable[j];
            ctable[j] = ctable[j+1];
            Ltable[j] = (ntable[j] * Ltable[j] + ntable[j+1] * Ltable[j+1])
                       / (ntable[j] + ntable[j+1]);
            ntable[j] += ntable[j+1];
            Evid += (ctable[j] - ctable[j-1]) * Ltable[j];
        }
        Common->Evidence = Evid;
        Common->Information = ctable[j] * Ltable[j] - Evid;

// Exit?
    } while( ! finish );
    CALLvalue = finish;

Exit:
// Free memory
    BayesFree(Common, Objects, Links, Ocean);
    FREE(Links)
    FREE(ntable)
    FREE(ctable)
    FREE(Ltable)
    return CALLvalue;
#undef NEXTRA
}

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function:  BayeShape
// 
// Purpose:   Transfer BayeSys hypercube to User's shape.
//
//   Shape  Description           Prob(Coord[i])                    Range of i
//
//     0   Permutation         uniform on integers Perm(0...N-1)     0...N-1
//     1   Positive orthant    exp(-x) in x > 0                      0...N-1
//     2   Simplex volume      uniform in SUM(x) < 1                 0...N-1
//     3   Simplex surface     uniform on SUM(x) = 1                  0...N
//     4   Ordered             uniform in 0 < x[0] < x[1] < ... < 1  0...N-1
//     5   Bell                Normal(0,1) in -inf < x < inf         0...N-1
//     6   Sphere volume       uniform in SUM(x^2) < 1               0...N-1
//     7   Hemisphere surface  uniform on SUM(x^2) = 1, x[N]>0        0...N
//     8   Sphere surface      uniform on SUM(x^2) = 1                0...N
//
// Notes:     The hypercube has 2^N corners, at which there may be singular
//            anisotropy in the exploration.  It is better to keep any such
//            singularities separated and diluted.
//
//     0   Permutation of integers 0,1,...,N-1.
//         No singularity, permutation given by ranked order of Cube[0...N-1].
//
//     1   Positive orthant
//         No singularity.  
//         The origin remains unchanged.
//         Other corners are at infinity.
//
//     2   Simplex volume
//         The origin remains unchanged.
//         Its N adjacent corners go to the remaining N vertices x[i]=1
//         of the simplex.
//         The remaining 2^N-N-1 corners are isolated singularities
//         distributed over the simplex roof SUM(x)=1.
//
//     3   Simplex surface
//         The origin moves to the upper vertex x[N]=1.
//         Its N adjacent corners go to the other surface vertices x[i]=1.
//         The remaining 2^N-N-1 corners are isolated singularities
//         distributed over the floor x[N]=0.
//
//     4   Ordered
//         The origin does not move.
//         Its N adjacent corners go to the other vertices (0,..0,0,1,1,..,1).
//         The remaining 2^N-N-1 corners are isolated singularities
//         distributed over the roof x[N-1]=1.
//         This is better than just sorting the x[i] from a hypercube,
//         with its N!-to-1 overlapping.
//
//     5   Bell
//         No singularity.
//         The hypercube centre becomes the origin.
//         All corners project out to infinity.
//
//     6   Sphere volume
//         All corners are isolated singularities on the sphere surface
//         SUM(x^2)=1, distributed with their original cubic symmetry.
//         This is much better than concentrating all the singularity at
//         the centre, as with using polar coordinates.
//
//     7   Hemisphere surface
//         All corners are isolated singularities at the equator x[N]=0,
//         distributed with their original cubic symmetry.
//
//     8   Sphere surface
//         All corners are isolated singularities at the surface's equator
//         x[N]=0, distributed with their original cubic symmetry.
//         The penalty for having this arrangement is that the north and south
//         hemispheres touch everywhere, with even Hilbert points defined as
//         north x[N]>0 and odd Hilbert points defined as south x[N]<0.
//         This will slow exploration because about half of the trial
//         points will be in the wrong hemisphere.  And there will be trouble
//         if an over-intelligent user tries to hide some information of his
//         own in this bit of lowest arithmetical significance.
//         The alternative, of collecting all the singularities together
//         down at the south pole, is worse.
//
// History:   JS        18 Sep 2003, 20 Oct 2003
//-----------------------------------------------------------------------------
//
void BayeShape(
double*    Coord,     //   O  coordinates              [N] or [N+1]
double*    Cube,      // I    hypercube position       [N]
int        N,         // I    dimension
int        Shape)     // I    choice of shape
{
static const double logG1[20] = { // log( (N/2)! )
            0.0000000000000000, -0.1207822376352452,  0.0000000000000000,
            0.2846828704729193,  0.6931471805599453,  1.2009736023470752,
            1.7917594692280550,  2.4537365708424432,  3.1780538303479458,
            3.9578139676187183,  4.7874917427820458,  5.6625620598571409,
            6.5792512120101012,  7.5343642367587300,  8.5251613610654147,
            9.5492672573009951, 10.6046029027452509, 11.6893334207972703,
           12.8018274800814691, 13.9406252194037652};
static const double  Z = (unsigned)(-1) + 1.0;          // 2^32
static const double  logpibytwo   =  0.45158270528945486473;

    unsigned  hemisphere;
    double    p, q, r, rr, s, t, inward, inwardN, outward, scale;
    int       i;

    switch( Shape )
    {
    case 0:
        ShapeIndex(N, Cube, Coord);
        break;

    case 1:
    case 2:
    case 3:
    case 4:
        for( i = 0; i < N; i++ )
            Coord[i] = -log(Cube[i]);
        if( Shape == 1 )
            break;

        r = 0.0;
        for( i = 0; i < N; i++ )
            r += Coord[i];                     // Old radius
        s = t = q = exp(-r/N);
        p = r * q;
        for( i = 1; i < N; i++ )
        {
            t *= p / i;
            s = s * q + t;
        }                                      // Outward cumulant
        s = (s * s < DBL_EPSILON) ? 1.0 - s / N
           : pow(1.0 - s, 1.0 / N);            // New radius
        s /= r;
        for( i = 0; i < N; i++ )               // Rescale
            Coord[i] *= s;
        if( Shape == 2 )
            break;

        t = 1.0;
        for( i = 0; i < N; i++ )
            t -= Coord[i];
        Coord[N] = t;
        if( Shape == 3 )
            break;

        for( i = 1; i <= N; i++ )
            Coord[i] += Coord[i-1];
        break;

    case 5:
    case 6:
    case 7:
    case 8:
        for( i = 0; i < N; i++ )
            Coord[i] = InvNorm(Cube[i]);
        if( Shape == 5 )
            break;

        rr = 0.0;
        for( i = 0; i < N; i++ )
        {
            Coord[i] = t = InvNorm(Cube[i]);
            rr += t * t;
        }
        r = sqrt(rr);
        scale = (N < 20) ? logG1[N] : logGamma(N/2. + 1.);
        scale = sqrt(2.0) / exp(scale / N);
        inwardN = r * scale;
        inward = pow(inwardN, N);
        outward = 1.0 - inward;                    // polar limit
        if( inward * inward > DBL_EPSILON )        // non-polar evaluation
        {
            if( N & 1 )
            {
                q = exp(-rr / (N + 1.0));
                p = q * rr;
                t = q / rr;
                s = 0.0;
                for( i = 1; i < N; i += 2 )
                {
                    t *= p / i;
                    s = s * q + t;
                }
                s = exp(logerf(r, 1) - 0.5 * (rr + logpibytwo)) + r * s;
            }
            else
            {
                t = s = q = exp(- rr / N);
                p = q * rr;
                for( i = 2; i < N; i += 2 )
                {
                    t *= p / i;
                    s = s * q + t;
                }
            }
            outward = s;
            inward = 1.0 - outward;
            inwardN = pow(inward, 1.0 / N);
        }
        if( Shape == 6 )
        {
            scale = (outward * outward < DBL_EPSILON)
                   ? (1.0 - outward / N) / r : inwardN / r;
            for( i = 0; i < N; i++ )
                Coord[i] *= scale;
            break;
        }

        ShapeSolve(N, outward, inwardN, &scale, &Coord[N]);
        scale /= r;
        for( i = 0; i < N; i++ )
            Coord[i] *= scale;
        if( Shape == 7 )
            break;

        hemisphere = 0;
        for( i = 0; i < N; i++ )
            hemisphere ^= (unsigned)(Cube[i] * Z);          // Hilbert parity
        if( hemisphere & 1 ) 
            Coord[N] = -Coord[N];
        break;

    default:
        for( i = 0; i < N; i++ )
            Coord[i] = Cube[i];
        break;
    }
}