Sampling Rosenbrock Density

This is some example code for sampling the Rosenbrock density following Goodman & Weare (2010).

The Rosenbrock density [ \ln , \pi (x_1,x_2) \propto -\frac{1}{20} \left [ 100(x_2 - x_1^2)^2 + (1-x_1)^2 \right ] ]

has the following contours:

Here's some code to sample the density.

import numpy as np
import pylab as pl

import markovpy as mp
import markovpy.diagnostics

The density function is:

def lnposterior(p):
    return -(100*(p[1]-p[0]*p[0])**2+(1-p[0])**2)/20.0

We start with the walkers sampled from a two dimensional Gaussian centered on ( (x_1,x_2) = (0,0) ) with variance ( 10^2 ) in both dimensions. This is a pretty bad guess but we'll find that that's no problem.

np.random.seed()
nwalkers = 100
pos0 = []
for k in range(nwalkers):
    pos0.append(10.*np.random.randn(2))
pos0=np.array(pos0)

Then, we grab the current random number state to pass to the sampler. In fact, this isn't necessary – the sampler will seed itself – but it shows how to pass a specific state if you require that.

state0 = np.random.get_state()

sampler = mp.ensemble.EnsembleSampler(nwalkers,2,lnposterior)

for position,prob,state in sampler.sample(pos0,None,state0,iterations=2000):
    pass

Then, the results can be plotted:

mp.plotting.plot2d(sampler.chain[:,0,:].flatten(),sampler2.chain[:,1,:].flatten(),bins=100)

It is interesting to note that an equivalent sampling would result from:

sampler2 = mp.ensemble.EnsembleSampler(nwalkers,2,lnposterior)
for position2,prob2,state2 in sampler2.sample(pos0,None,state0,iterations=500):
    print_some_stuff()

run_some_diagnostics_or_something_if_you_want()

for position2,prob2,state2 in sampler2.sample(position2,prob2,state2,iterations=1000):
    print_some_other_stuff()

but it allows manipulations of the sample at any intermediate step.