Hi,
I want to implement a new sampler and need the banana-shaped distribution as a test case.
Here I found an implementation for pymc, but I can’t transfer it to pymc3.
This is how I can transform samples from a multivariate normal into the banana-shape:
def banana(X, b=0.03):
"""Twist the second column of X into a banana."""
X = [x for x in X.copy()]
X[1] += b * X[0]**2 - 100 * b
return X
C = [
[100,1],
[1,1]
]
bananasamples = banana(numpy.random.multivariate_normal([0, 0], C, size=200).T)
Then it comes to the pymc3 model, I have a fundamental misunderstanding, but I can’t figure out what it is.
def run_stepper(stepper_cls):
with pymc3.Model() as pmodel:
x = pymc3.Uniform('x', lower=-50, upper=50)
y = pymc3.Uniform('y', lower=-50, upper=50)
G = pymc3.MvNormal('G', mu=[0,0], cov=C)
# ???????
pymc3.Potential('X', G.logp([x,y]))
step = stepper_cls()
mt = pymc3.sample(step=step, tune=1000)
return mt
if __name__ == '__main__':
methods = [pymc3.Metropolis, pymc3.NUTS]
fig, hosts = visualization.fig_hosts(1, len(methods), 5, 5)
for c in range(len(methods)):
mt = run_stepper(methods[c])
hosts[0,c].scatter(S[0], S[1], s=0.2)
hosts[0,c].scatter(mt['x'], mt['y'], s=0.3)
accept_stat = [s for s in mt.stat_names if 'accept' in s][0]
accept_rate = mt.get_sampler_stats(accept_stat).mean()
title = '{} ({:.1f}%)'.format(methods[c].name, accept_rate * 100)
hosts[0,c].set_title(title)
pyplot.show()
I am sure there is a one- or two-line way to do it right and would appreciate any help.
The expectation then is that Metropolis will of course have a much worse acceptance rate compared to NUTS.
thanks
michael
