2013-07-11 4 views
7

이것은 아마도 어리석은 질문 일 것입니다.사용자 정의 PyMC 배포 정의

저는 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 또는 이것 저것 같은 것을 사용하지만, 나는 주석 밀도 함수에 적합 할 찾은 예.

도움을 주시면 감사하겠습니다.

답변

9

쉬운 방법은 확률 장식 사용하는 것입니다 내 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을 사용하여 커널 밀도 추정기가있는 노드를 만드는 예제가 포함되어 있습니다.

+0

대단히 유용했습니다. 감사합니다. – stellographer

+0

@jcrudy. 내 자신의 간단한 사전 정의하려고 할 때 귀하의 대답을 실행합니다. 이 질문을 피하기 위해, 나는 내 자신의 [여기] (http://stackoverflow.com/questions/23198247/custom-priors-in-pymc)를 시작했고, 당신이 그것에 빛을 비추는 지 궁금했다. 감사. –

관련 문제