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.