Defining a custom PyMC distribution

Question

Perhaps this is a stupid question.

I am trying to fine-tune data into a very strange PDF file using the MCMC rating in PyMC. In this example, I just want to figure out how to fit the normal distribution, where I manually enter a regular PDF. My code is:

data = []; for count in range(1000): data.append(random.gauss(-200,15)); mean = mc.Uniform('mean', lower=min(data), upper=max(data)) std_dev = mc.Uniform('std_dev', lower=0, upper=50) # @mc.potential # def density(x = data, mu = mean, sigma = std_dev): # return (1./(sigma*np.sqrt(2*np.pi))*np.exp(-((x-mu)**2/(2*sigma**2)))) mc.Normal('process', mu=mean, tau=1./std_dev**2, value=data, observed=True) model = mc.MCMC([mean,std_dev]) model.sample(iter=5000) print "!" print(model.stats()['mean']['mean']) print(model.stats()['std_dev']['mean']) 

In the examples I found, everyone uses something like mc.Normal or mc.Poisson or something else, but I want to match a function with density comments.

Any help would be appreciated.

Answer 1

A simple way is to use a stochastic decorator:

 import pymc as mc import numpy as np data = np.random.normal(-200,15,size=1000) mean = mc.Uniform('mean', lower=min(data), upper=max(data)) std_dev = mc.Uniform('std_dev', lower=0, upper=50) @mc.stochastic(observed=True) def custom_stochastic(value=data, mean=mean, std_dev=std_dev): return np.sum(-np.log(std_dev) - 0.5*np.log(2) - 0.5*np.log(np.pi) - (value-mean)**2 / (2*(std_dev**2))) model = mc.MCMC([mean,std_dev,custom_stochastic]) model.sample(iter=5000) print "!" print(model.stats()['mean']['mean']) print(model.stats()['std_dev']['mean']) 

Note that my custom_stochastic function returns the probability of registration, not probability, and that this is the log probability for the entire sample.

There are several other ways to create custom stochastic nodes. This doc gives more details, and this gist contains an example using pymc.Stochastic to create a node with kernel density estimation.