カスタムPyMCディストリビューションを定義する

Question

これはおそらく愚かな質問です。

私は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やその他もろもろのようなものを使用しますが、私はコメントアウト密度関数に収まるようにしたい見つけた例。

ご協力いただければ幸いです。

Answer 1

簡単な方法は、確率的なデコレータを使用することです：私のcustom_stochastic関数は対数尤度ではなく、可能性を返す

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

注意をし、それが全体のサンプルに対する対数尤度であること。

カスタム確率的ノードを作成するには、他にもいくつかの方法があります。このdocは詳細を示し、gistには、pymc.Stochasticを使用してカーネル密度推定子を持つノードを作成する例が含まれています。