DNest (the Diffusive Nested Sampler) requires the GNU scientific 
library. To compile on a GNU/Linux system, simply run 'make'. You'll 
need g++ of course. This will make an executable called main, and a 
shared library in case you prefer to link to DNest from your own code 
rather than just copying and modifying the DNest code.

To implement your own problem, edit TestModel.h and TestModel.cpp to 
define the parameters, Metropolis-Hastings proposal distribution (edit 
TestModel::perturb()), prior distribution (currently considered implicit 
in the proposal distribution, but this can be modified) and likelihood 
function. Numerical parameters for the algorithm live in OPTIONS and are 
loaded at runtime. Sensible defaults have been provided. The current 
version has no termination condition. It will run forever - make sure 
you kill it when you have enough output! To view the results, run 
showresults.py using Python2. You will need numpy and pylab. This script 
computes the posterior weights, a posterior sample, and the evidence.

sample.txt contains all of the samples, with one per line. However, 
these samples are not from the posterior. To make the samples represent 
the posterior, you need to weight them appropriately. The showresults 
script calculates these weights, and outputs resampled (equally 
weighted) samples to posterior_sample.txt.

> The next step will be to figure out the right way to code up my 
> likelihood function under "void TestModel::calculateLogLikelihood()" 
> in TestModel.cpp ...

This shouldn't be too difficult. Firstly, within TestModel there is a 
vector (array) called x, which is currently 20 elements. You should find 
the constructor and change the number 20 to however many parameters you 
have.

TestModel::TestModel()
:x(20) // <---------change this number
{

}

Then, you need to modify void TestModel::fromPrior() such that the 
computer can generate random examples of your model from the prior 
distribution. For example the current fromPrior() looks like this:

void TestModel::fromPrior()
{
	for(size_t i=0; i<x.size(); i++)
		x[i] = -0.5 + randomU();
	Model::fromPrior();
	calculateLogLikelihood();
}

The first two lines generate all of the x's independently from a uniform 
distribution between -0.5 and 0.5, this is what you'll need to change. 
Use the functions randomU() and randn() to generate from uniform(0,1) 
and standard normal(0, 1) respectively (these are not built in C++ 
functions, I define them in RandomNumberGenerator which calls gsl). The 
last two lines *must* be retained.

Next thing you need to modify is TestModel::perturb(). This defines 
proposal distributions for Metropolis-Hastings exploration. At the 
moment I choose one of the parameters at random and perturb it by a 
randn() whose standard deviation is between 10^(-5) and 10^(1):

x[which] += pow(10.0, 1-6*randomU())*randn();

This randomised step size is a good idea, it's more efficient than Slice 
Sampling! . You may need to tune the range of stepsizes for your 
problem. The next line (involving mod) just wraps the proposed value 
into [-0.5 0.5] which was the domain for the test problem. You probably 
won't want to do that. The return value for this function accounts for 
the case when the proposal does not implicitly explore the prior for 
your problem, in which case you need to return log(prior_new) - 
log(prior_old). For example, if you have random walk proposals that 
explore a uniform distribution, but you actually want to explore a 
normal distribution.

You need to retain the lines

Model::perturb();
calculateLogLikelihood();

in TestModel::perturb().

Finally, you need to implement calculateLogLikelihood(). This can do 
whatever it needs to do, as long as the variable logl gets assigned at 
the end, you're good to go. If likelihood=0 anywhere, set logl=-1E300 
because C++ doesn't handle -infinity well at all.

I hope this helps - good luck with it!

-- Dr Brendon J. Brewer
Department of Physics, University of California, Santa Barbara
Web: http://www.physics.ucsb.edu/~brewer/

