gibbs1d.h:
#ifndef __GIBBS1D_H__
#define __GIBBS1D_H__
#include <vector>
/**
* Produces a sequence of variates with P(x) propto exp(-beta * H(nodes))
* as its stationary distribution, where H(nodes) = the number of consecutive
* pairs with equal state.
*/
class Gibbs1d {
public:
/**
* beta -- larger values imply larger affinities for forming clusters.
* numNodes -- number of points along the line.
* q -- number of states for each point.
*/
Gibbs1d(double beta, int numNodes, int q = 2);
/**
* Update the state of each node by Gibbs sampling.
*/
void iterate();
/**
* Returns a pointer to a numNodes length array giving
* the current state of each node.
*/
const int* current() const;
private:
//Hide the copy ctor and assignment operator
Gibbs1d(const Gibbs1d&);
Gibbs1d& operator=(const Gibbs1d&);
const double beta;
const int numNodes;
const int q;
std::vector<int> states;
const double probNotEq1;
const double probNotEq2;
const double probEq;
const double prob1Out;
};
#endif
gibbs1d.cpp:I think there must be a much more elegant implementation of this. For one thing, the use of "sentinel" values isn't providing much of a payoff here. And the code for computing the probabilities is just very ugly. What's a better alternative to the nested if-statements?
#include <algorithm>
#include <cmath>
#include "gibbs1d.h"
#include "mtrandom.h"
using namespace std;
/**
* beta -- larger values imply larger affinities for forming clusters.
* numNodes -- number of points along the line.
* q -- number of states for each point.
*/
Gibbs1d::Gibbs1d(double beta, int numNodes, int q)
: beta(beta),
numNodes(numNodes),
q(q),
states(numNodes + 2),
probNotEq1(exp(beta) / (2.0 * exp(beta) + q - 2)),
probNotEq2(2.0 * probNotEq1),
probEq(exp(2.0 * beta) / (exp(2.0 * beta) + q - 1)),
prob1Out(exp(beta) / (exp(beta) + q - 1))
{
//Add sentinel values to the ends to avoid special casing
//the boundaries.
states[0] = states[numNodes + 1] = -1;
}
/**
* Returns an integer in the interval [0, q-1]
* != exclude, assigning equal probability to each
* allowed value.
* Assumes exclude is in [0, q-1].
*/
static int uniformNotEq(int q, int exclude)
{
int r = (int) ((q - 1) * genrand_real2());
if (r >= exclude) {
return r + 1;
}
else {
return r;
}
}
/**
* Returns an integer in the interval [0, q-1]
* != exclude1 or exclude2, assigning equal probability to each
* allowed value.
* Assumes exclude1 and exclude2 are in [0, q-1]
* and that exclude1 != exclude2.
*/
static int uniformNotEq(int q, int exclude1, int exclude2)
{
if (exclude1 > exclude2) {
swap(exclude1, exclude2);
}
int r = (int) ((q - 2) * genrand_real2());
if (r >= exclude1) {
r = r + 1;
}
if (r >= exclude2) {
r = r + 1;
}
return r;
}
/**
* Returns the given state with the given probability.
* Otherwise, chooses one of the other q - 1 states
* uniformly.
*/
static int withProbElseUniformNotEq(double prob, int state, int q)
{
const double u = genrand_real2();
if (u <= prob) {
return state;
}
else {
return uniformNotEq(q, state);
}
}
void Gibbs1d::iterate()
{
for (int i = 1; i <= numNodes; ++i) {
const int left = states[i - 1];
const int right = states[i + 1];
if (left >= 0) {
if (right >= 0) {
if (left != right) {
const double u = genrand_real2();
if (u <= probNotEq1) {
states[i] = left;
}
else if (u <= probNotEq2) {
states[i] = right;
}
else {
states[i] = uniformNotEq(q, left, right);
}
}
else { // left == right
states[i] = withProbElseUniformNotEq(probEq, left, q);
}
}
else { //right < 0
states[i] = withProbElseUniformNotEq(prob1Out, left, q);
}
}
else {
if (right >= 0) {
states[i] = withProbElseUniformNotEq(prob1Out, right, q);
}
else {
states[i] = (int) (q * genrand_real2());
}
}
}
}
/**
* Returns a pointer to a numNodes length array giving
* the current state of each node.
*/
const int* Gibbs1d::current() const
{
return &states[1];
}
//I really hate how Blogger cuts off the edges of wide elements. I suppose I can edit the style to fix this.
