gibbssampler.cpp

gibbs

#include "GibbsSampler.h"
#include "VarConfig.h"
#include <algorithm>
#include "Prob.h"

#define CHECKMEM(x) \
{ \
    cout << "Pre-" << x << endl; \
    int* foo = new int[1000]; \
    cout << "Post-" << x << endl; \
}

#define SQR(x) ((x) * (x))

#if 0
// Unused -- we use the VarConfig class instead.
int GibbsSampler::configIndex(VarSet& allVars, list<int>& testIndices)
{
    int rangeProduct = 1;
    int ret = 0;
    list<int>::iterator i;
    for (i = testIndices.begin(); i != testIndices.end(); i++) {
        // Test vars must be discrete.
        ret += rangeProduct * (int)allVars[*i];
        rangeProduct *= model.getRange(*i);
    }

    return ret;
}
#endif


// Old way of burning in -- run a constant number of iterations.
void GibbsSampler::burnInChain(VarSet &chain, const VarSet& evidence, 
        int burnInIters) const
{
    for (int iter = 0; iter < burnInIters; iter++) {
        for (int v = 0; v < evidence.getNumVars(); v++) {
            if (!evidence.isTested(v)) {
                chain[v] = model.MBsample(v, chain);
            }
        }
    }
}

unsigned int GibbsSampler::sampleFromDist(const vector<double>& dist) const
{
    double p = (double)rand()/RAND_MAX;
    for (unsigned int v = 0; v < dist.size(); v++) {
        p -= dist[v];
        if (p < 0.0) {
            return v;
        }
    }

    // DEBUG
    cout << "Error: returning last value in sampleFromDist()\n";

    return (dist.size() - 1);
}

double GibbsSampler::testConvergence(vector<vector<double> > summaries,
        vector<vector<double> > sqSummaries, int n) const
{
    double maxR = 0.0;
    int numChains = summaries.size();
    int numSummaries = sqSummaries[0].size();
    vector<double> avgSummaryVariance(numSummaries);

    // Iterate through all chains
    for (int c = 0; c < numChains; c++) {
        // Compute within-chain variance for each summary statistic
        for (int s = 0; s < numSummaries; s++) {

#if 0
            // HACK: We require that every single state receive at least a
            // fractional count in some chain before we converge.
            if (summaries[c][s] == 1.0) {
                // DEBUG
                cout << "c = " << c << "; s = " << s << endl;
                return 1000.0;
            }
#endif

            double chainVar = (sqSummaries[c][s] - SQR(summaries[c][s])/n)/(n-1);
            avgSummaryVariance[s] += chainVar/numChains;
        }
    }

    for (int s = 0; s < numSummaries; s++) {

        // Compute between-chain variance for each summary statistic
        double squareSum = 0.0;
        double sum = 0.0;

        for (int c = 0; c < numChains; c++) {
            squareSum += SQR(summaries[c][s]/n);
            sum += summaries[c][s]/n;
        }
        double betweenChainVariance
            = (squareSum - SQR(sum)/numChains)/(numChains-1);

        // Compute convergence criteria R
        double R = ((n-1.0)/n*avgSummaryVariance[s] 
                     + betweenChainVariance)/avgSummaryVariance[s];

        // Report largest convergence statistic
        if (R > maxR) {
            maxR = R;
#if 0
        // DEBUG
        cout << "Counts:";
        for (int c = 0; c < numChains; c++) {
            cout << " " << summaries[c][s];
            cout << " (" << sqSummaries[c][s] << ")";
        }
        cout << "\n";
        cout << "Within: " << avgSummaryVariance[s];
        cout << "; Between: " << betweenChainVariance;
        cout << "; R: " << sqrt(R) << endl;
        // END DEBUG
#endif
        }
    }

    // DEBUG
    //cout << "sqrt(R) = " << sqrt(maxR) << endl;

    return sqrt(maxR);
}

double GibbsSampler::predictIters(vector<vector<double> > summaries,
        vector<vector<double> > sqSummaries, int n) const
{
    double maxV = 0.0;
    int numChains = summaries.size();
    int numSummaries = sqSummaries[0].size();
    vector<double> avgSummaryVariance(numSummaries);

    // Iterate through all summary statistics
    for (int s = 0; s < numSummaries; s++) {

        // Consider the chains as independent estimates of each summary
        // statistic, and compute their standard deviation.
        double squareSum = 0.0;
        double sum = 0.0;

        for (int c = 0; c < numChains; c++) {
            squareSum += SQR(summaries[c][s]/n);
            sum += summaries[c][s]/n;
        }

        // See page 740 of DeGroot and Schervish, 3rd ed.
        double S = sqrt((squareSum - SQR(sum)/numChains)/numChains);
        double sigma_hat = sqrt((double)n) * S;

        // Compute number of expected iterations (with 95% certainty)
        // to get the estimate correct within 5%.
        // (See page 707, eqn 11.1.5 of DeGroot and Schervish, 3rd ed.)
        double epsilon = 0.05 * sum/numChains;
        double v = SQR(1.96 * sigma_hat/epsilon);

        // Report largest number of iterations to run
        if (v > maxV) {
            maxV = v;

#if 0
            // DEBUG
            cout << "epsilon = " << epsilon << endl;
            cout << "sigma_hat = " << sigma_hat << endl;
            cout << "n = " << n << endl;
            cout << "sum = " << sum << endl;
            cout << "squareSum = " << squareSum << endl;
#endif

#if 0
            double mean = sum/numChains;
            double maxRatio = 1.0;
            for (int c = 0; c < numChains; c++) {
                double currStat = summaries[c][s]/n;
                if (currStat/mean > maxRatio) {
                    maxRatio = currStat/mean;
                }
                if (mean/currStat > maxRatio) {
                    maxRatio = mean/currStat;
                }
            }
            cout << "Max ratio: " << maxRatio << endl;
#endif
        }
    }

    return maxV;
}


void GibbsSampler::runMarginalInference(const VarSet& evidence)
{
    // We use this vector to convert var/value pairs into summary 
    // statistic indices.
    vector<vector<int> > index(model.getNumVars());
    int numSummaries = 0;
    for (int v = 0; v < model.getNumVars(); v++) {
        for (int val = 0; val < model.getRange(v); val++) {
            index[v].push_back(numSummaries++);
        }
    }

    // Keep track of marginal counts for all test variables
    vector<vector<double> > counts(numChains, vector<double>(numSummaries));
    vector<vector<double> > sqCounts(numChains, vector<double>(numSummaries));
    for (int c = 0; c < numChains; c++) {
        for (int s = 0; s < numSummaries; s++) {
            counts[c][s] = 0.0;
            sqCounts[c][s] = 0.0;
            //counts[c][s] = 1.0;
            //sqCounts[c][s] = 1.0;
        }
    }

    // Initialize and burn-in all chains
    vector<VarSet> chains(numChains);
    for (int c = 0; c < numChains; c++) {
        chains[c] = evidence;
        model.wholeSample(chains[c]);

        // Use a fixed number of burn-in iters, if appropriate
        if (fixedIters) {
            burnInChain(chains[c], evidence, burnInIters);
        }
    }

    // Sample, sample, sample until convergence
    double burnin_iter = 0;
    double sampling_iter = 0;
    double predicted_iters = minIters;
    bool burnin_done = fixedIters;
    while (1) {
        if (burnin_done) {
            sampling_iter++;
        } else {
            burnin_iter++;
        }

        // Sample all variables and increase counts
        for (int c = 0; c < numChains; c++) {
            for (int v = 0; v < model.getNumVars(); v++) {

                // Don't resample evidence variables
                if (evidence.isTested(v)) {
                    continue;
                }
                
                // Sample
                vector<double> dist = model.MBdist(v, chains[c]);
                chains[c][v] = sampleFromDist(dist);

#define RAOBLACKWELL
#ifdef RAOBLACKWELL
                // Update (Rao-Blackwellized) counts
                for (int val = 0; val < model.getRange(v); val++) {

                    // Update counts using the distribution
                    int s = index[v][val];
                    counts[c][s] += dist[val];
                    sqCounts[c][s] += dist[val]*dist[val];
                }
#else
                // Update counts
                int s = index[v][(int)chains[c][v]];
                // DEBUG
                //cout << "s = " << s << endl;
                counts[c][s]++;
                sqCounts[c][s]++;
#endif
            }
        }

        // After completing some minimum number of iterations,
        // check for convergence of our burn-in period
        if (!burnin_done && burnin_iter >= burnInIters 
                && ((int)burnin_iter % 100 == 0) 
                && (testConvergence(counts, sqCounts, (int)burnin_iter)
                    < convergenceRatio)) {

            // Stop burn-in
            burnin_done = true;

            // Throw away counts for burn-in period
            for (int c = 0; c < numChains; c++) {
                for (int s = 0; s < numSummaries; s++) {
                    counts[c][s] = 0.0;
                    sqCounts[c][s] = 0.0;
                    //counts[c][s] = 1.0;
                    //sqCounts[c][s] = 1.0;
                }
            }
        }


        // Test for convergence of the sampling
        // Go until our standard error among the different chains is
        // less than 5% of the predicted value.
        if (burnin_done && sampling_iter >= predicted_iters) {

            // Stop, if we're only running a fixed number of iterations
            if (fixedIters) {
                break;
            }
            
            predicted_iters = predictIters(counts, sqCounts, (int)sampling_iter);

            // DEBUG
            cout << "Predicted iters = " << predicted_iters << endl;

            if (predicted_iters <= sampling_iter) {
                break;
            }
        }
    }

    // Save distributions
    for (int v = 0; v < model.getNumVars(); v++) {

        if (evidence.isTested(v)) {
            continue;
        }
        Distribution m(model.getRange(v));
        for (int val = 0; val < model.getRange(v); val++) {
            m[val] = 0;
            for (int c = 0; c < numChains; c++) {
                m[val] += counts[c][index[v][val]];
            }
        }
        m.normalize();
        marginals[v] = m;

#ifdef DEBUG
        if (!evidence.isTested(v)) {
            cout << v << ": " << m << endl;
        }
#endif
    }

    // DEBUG
    if (!fixedIters) {
        cout << burnin_iter << "; " << sampling_iter << endl;
    }
}

void GibbsSampler::runJointInference(const list<int>& queryVars,
        const VarSet& evidence)
{
    VarSchema schema = model.getSchema();
    VarConfig query(evidence, queryVars, schema);
    int numSummaries = query.getMaxIndex() + 1;
    
    // Keep track of counts for all test configurations
    vector<vector<double> > counts(numChains, vector<double>(numSummaries));
    vector<vector<double> > sqCounts(numChains, vector<double>(numSummaries));
    for (int c = 0; c < numChains; c++) {
        for (int s = 0; s < numSummaries; s++) {
            counts[c][s] = 1.0/numSummaries;
            sqCounts[c][s] = 1.0/numSummaries;