import numpy as np
import emcee

def lnprob(x,mu,icov):
    diff = x - mu
    return -np.dot(diff,np.dot(icov,diff))/2.0

ndim = 50

means = np.random.rand(ndim)

cov = 0.5 - np.random.rand(ndim ** 2).reshape((ndim,ndim))
cov = np.triu(cov)
cov += cov.T - np.diag(cov.diagonal())
cov = np.dot(cov,cov)

icov = np.linalg.inv(cov)

nwalkers = 250
p0 = np.random.rand(ndim * nwalkers).reshape((nwalkers,ndim))

sampler = emcee.EnsembleSampler(nwalkers,ndim,lnprob,arg=[means,icov])


lnprob(p,means,icov)
pos, prob, state = smapler.run_mcmc(p0,100)
sampler.reseet()

sampler.run_mcmc(pos,1000)


import matplotlib.pyplot as pl

for i in range(ndim):
    pl.figure()
    pl.hist(sampler.flatchain[:i],100, color="k", histtype = "step")
    pl.title("Dimension {0,d}".format(i))

pl.show()

print("Mean acceptance fraction:  {0:.3f}").format(np.mean(sampler.acceptance_fraction))
