# PyMC: how can I use a custom sampler for one specific variable in a model?

Basically I have a function to efficiently sample one variable in my model that I would like PyMC to use, along with its normal samplers for everything else. Since I have this function I don't have a need to define a log likelihood for the variable either.

Can I do this using the Stochastic class and making a random function? Or do I need to make a custom step method?

I think you will need to make a custom step method. This is what I did for a project on sampling random spanning trees with PyMC2, and you can see my spanning tree Metropolis step method here.

Here is a minimal PyMC2 version:

class StandardNormal(pm.Gibbs):
def __init__(self, stochastic, verbose=None):
pm.Gibbs.__init__(self, stochastic, verbose=verbose)

def step(self):
self.stochastic.value = np.random.normal()

A = pm.Uninformative('A', value=0)
B = pm.Normal('B', mu=20+A, tau=10)

mc = pm.MCMC([A,B])
mc.use_step_method(StandardNormal, A)
mc.sample(iter=5000, burn=100)


And here is a corresponding PyMC3 version:

class StandardNormal(object):
def __init__(self, var):
self.var = var.name

def step(self, point):
new = point.copy()
new[self.var] = np.random.normal()

return new

with pm.Model() as model:
A = pm.Flat('A')
B = pm.Normal('B', mu=20+A, sd=10)

step_A = StandardNormal(var=A)
step_B = pm.step_methods.Metropolis(vars=[B])

trace = pm.sample(5000, [step_A, step_B])


You can see it all in action here.

• Thanks, I also just came across this (written by you I believe!) which seems relevant and fairly straightforward. Any idea if it is compatible with PyMC3? I'm still a bit surprised this isn't as simple as providing a Stochastic with just a random function and having PyMC recognize that automatically. May 9 '14 at 16:42
• Followup question: if I used a custom step method for sampling one of my variables, then how do I even define the variable in PyMC given that I don't want to bother with a log likelihood function (since I'm defining my own sampler)? May 9 '14 at 18:07
• It looks like you got the answer already on github, but I will update here in case future searchers find this page. May 10 '14 at 0:14