1d Gauss example -- returns log(Z)=-0.32, instead of 0?

[JohannesBuchner]

Dear Brendon,

I would like to use your Diffusive Nested Sampling technique. To start, I created an extremely simple example -- a 1d Gaussian. However, I get incorrect results: It returns log(Z)=-0.32, instead of 0. Could you please help me figure out what I did wrong? Also, do you have advice on how to build the perturb() function for a given problem? So far, I reused the multi-scale loguniform/Gaussian from SpikeSlab. I am also unsure in which cases not to return 0.

Kind regards, Johannes

[eggplantbren]

Reflecting off the walls could be cool. Wrapping usually works fine though. The code in my answer could be dropped in for that.

Another issue might be the 'number of levels' in the OPTIONS file. For such a simple problem you only need about 5 compressions to get from the prior to the posterior. If you try to use a lot more, strange issues might occur due to havig a very narrow distribution which rejects all of your moves.

[eggplantbren]

The return value of 'perturb' is used to implement non-uniform priors. For example, here is an example with a Normal(0, 1) prior. The move I use implies a uniform prior, but you can compensate for that by returning the log of something that will get included in the acceptance probability.

void MyModel::fromPrior()
{
    x = randn();
}

double MyModel::perturb()
{
    double logH = 0.;

    logH -= -0.5*x*x;
    x += pow(10., 1.5 - 6*randomU())*randn();
    logH += -0.5*x*x;

    return logH;
}

[JohannesBuchner]

Yes, indeed, it was the OPTIONS. I had to reduce both the number of levels and the force. I am using this now:

1       # Number of particles
1000    # new level interval
100000  # save interval
200     # threadSteps - how many steps each thread should do independently before communication
5       # maximum number of levels
10      # Backtracking scale length (lambda in the paper)
5       # Strength of effect to force histogram to equal push. 0-10 is best. (beta in the paper)
200     # Maximum number of saves (0 = infinite)

I intend to do a comparison of various algorithms. Do you have some general advise on how to choose the parameters for a fair comparison?

For example for these problems:

• Gaussian in 10d (relatively small space)
• 2 thin rings in 2d (2 modes, very thin regions)
• Eggbox problem (~20 equal modes)

To this end, I started a Python interface to your software, which you can find at https://github.com/JohannesBuchner/PyDNest .

Cheers, Johannes

[JohannesBuchner]

By the way, if I interpret the FitSine example correctly, Diffusive Nested Sampling can be used to determine the dimension of models simultaneously with their parameters, similar to what RJMCMC should do. Or is this general to Nested sampling?

[eggplantbren]

Hi Johannes,

The gaussian in 10D might make Diffusive Nested Sampling look bad, because it's one where alternative methods work really well. ;-) In low dimensions you won't need so many levels. Just check there is a peak in the posterior weights plot.

The FitSine example is a trans-dimensional example. What it does is reversible jump MCMC. Diffusive Nested Sampling just asserts a different target distribution - not the posterior but a mixture of constrained priors.

• Brendon

[eggplantbren]

Thanks for PyDNest. That looks very useful! Many people would be confident implementing their model in Python but not C++. I wonder if it will work with the multithreading, though.

[JohannesBuchner]

Thank you for taking the time to explain things to me.

[JohannesBuchner]

I added a python re-implementation to PyDNest. The output should be compatible, but I am not convinced it is bug-free at this point.

[eggplantbren]

I am impressed! The ideas aren't that complicated, but whenever I've tried to implement it in another language I gave up because it became too complex. I'll let you know if I notice any issues.