Example using custom step method on PyMC: - generate samples for a sum of gaussian variables with mean 0 and std 1. - build a pymc model to fit a generative model to this set of observations. First declare all variables within the model context. The std of the (s variable) should be 1.4 as expected. Second try to update the variable p using a cu…

test_step_method.py

import pymc as mc
import numpy as np
import time
import theano.tensor as T
import theano
theano.config.compute_test_value = 'off'


## generate some data
N = 10000
obs_z = np.random.normal(size=N).reshape(N/2,2)
# this is a sum of two gaussians N(0,1) so it should be a N(0,1.4)
obs_z_sum = np.sum(obs_z,1)
print 'Observed', obs_z_sum.mean(), obs_z_sum.std()

#### Use a model with PyMC to fit to this data
with mc.Model() as model:
	z = mc.Normal('z',mu=0,sd=1,shape=(1,2)) 
	s = mc.Uniform('s',0,10)
	p = T.sum(z)
	z_sum = mc.Normal('z_sum',mu=p,sd=s,observed=obs_z_sum)
	step = mc.Metropolis() #mc.NUTS()
	traces = mc.sample(500, step,progressbar=True)

# will only print summary for s, as expected with mean 1.4
print mc.summary(traces,start=300)


#### Fit model using a custom step method
class MyExactUpdateStepMethod(object):
	    def step(self, point): #point is dict of parameter values
			 newpoint = point.copy()
			 newpoint['p'] = T.sum(newpoint['z'])
			 #update values in newpoint to be the exact formulas you talked about.
			 return newpoint

with mc.Model() as model:
	z = mc.Normal('z',mu=0,sd=1,shape=(1,2)) 
	s = mc.Uniform('s',0,10)
	p = 0.0 #T.sum(z)
	z_sum = mc.Normal('z_sum',mu=p,sd=s,observed=obs_z_sum)
	nuts_step = mc.Metropolis()
	exact_step = MyExactUpdateStepMethod()
	traces = mc.sample(500, step = [nuts_step,exact_step] ,progressbar=True)

print mc.summary(traces,start=300)