这可能是一个愚蠢的问题。

我试图使用PyMC中的MCMC评估来拟合数据到一个非常奇怪的PDF。对于这个例子，我只想弄清楚如何拟合到一个我手动输入的正态PDF。我的代码是：

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'])

我找到的例子都使用了类似mc.Normal或mc.Poisson之类的东西，但我想要拟合到注释掉的密度函数。

任何帮助都将不胜感激。

回答：

一种简单的方法是使用stochastic装饰器：

import pymc as mcimport numpy as npdata = 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'])

请注意，我的custom_stochastic函数返回的是对数似然，而不是似然，并且它是对整个样本的对数似然。

还有其他几种创建自定义随机节点的方法。这份文档提供了更多细节，而这个代码片段包含了一个使用pymc.Stochastic创建带有核密度估计的节点的例子。