In [1]:
import numpy as np
import matplotlib.pyplot as pl
import seaborn as sn
import emcee
import pymultinest
%matplotlib inline
In [2]:
class productTwoNumbers(object):
def __init__(self, y, noise):
self.y = y
self.noise = noise
self.lower = np.asarray([0.0, 0.0])
self.upper = np.asarray([50.0, 1.0])
def Prior(self, cube, ndim, nparams):
for i in range(ndim):
cube[i] = cube[i] * (self.upper[i]-self.lower[i]) + self.lower[i]
def Loglike(self, cube, ndim, nparams, lnew):
model = cube[0] * cube[1]
logL = -(self.y-model)**2 / (2.0*self.noise**2)
return logL
def logPosterior(self, theta):
if ( (np.any(theta < self.lower)) | (np.any(theta > self.upper)) ):
return -np.inf
else:
model = theta[0] * theta[1]
logL = -(self.y-model)**2 / (2.0*self.noise**2)
return logL
def sampleMultinest(self):
pymultinest.run(self.Loglike, self.Prior, 2, resume = False, verbose = False)
def sample(self):
ndim, nwalkers = 2, 100
self.theta0 = np.asarray([self.y/0.5,0.5])
#self.p0 = []
#for i in range(nwalkers):
#new = np.random.rand(2)*(self.upper-self.lower) + self.lower
#self.p0.append(new)
self.p0 = [self.theta0 + 0.01*np.random.randn(ndim) for i in range(nwalkers)]
self.sampler = emcee.EnsembleSampler(nwalkers, ndim, self.logPosterior)
self.sampler.run_mcmc(self.p0, 5000)
result = 10.0
noise = 1.0
out = productTwoNumbers(result, noise)
In [ ]:
out.sampleMultinest()
In [168]:
f, ax = pl.subplots(ncols=2, nrows=2, figsize=(12,8))
ax[0,0].plot(aSamples)
ax[0,1].hist(aSamples)
ax[0,0].set_ylabel('a')
ax[0,1].set_xlabel('a')
ax[1,0].plot(bSamples)
ax[1,1].hist(bSamples)
ax[1,0].set_ylabel('b')
ax[1,1].set_xlabel('b')
pl.tight_layout()
In [169]:
sn.set_style("dark")
aX = np.linspace(0, 50, 100)
bX = np.linspace(0, 1.0, 100)
X,Y = np.meshgrid(aX, bX)
Z = np.exp(-(result-X*Y)**2 / (2.0*noise**2))
f, ax = pl.subplots(figsize=(8,8))
ax.pcolor(X, Y, Z)
Out[169]:
In [170]:
pl.hexbin(aSamples,bSamples)
Out[170]:
In [76]:
p0 = [np.zeros(2) + 1*np.random.randn(2) for i in range(4)]
In [77]:
p0
Out[77]:
In [ ]: