pymc3 : 다중 관측 값

Question

필자는 매개 변수를 추정하고자하는 관찰 데이터를 가지고 있으며 PYMC3을 시험해 볼 수있는 좋은 기회라고 생각했습니다.

내 데이터는 일련의 레코드로 구성됩니다. 각 레코드는 고정 된 1 시간 기간과 관련된 한 쌍의 관찰을 포함합니다. 하나의 관찰은 주어진 시간 동안 발생하는 총 이벤트 수입니다. 다른 관찰은 그 기간 내의 성공 횟수입니다. 예를 들어, 데이터 포인트는 주어진 1 시간 동안 총 1000 개의 이벤트가 있었고 1000, 100 개의 이벤트는 성공한 이벤트라고 지정할 수 있습니다. 다른 기간에는 총 1000000 개의 이벤트가있을 수 있으며 그 중 120000 개가 성공한 이벤트입니다. 관측 값의 분산은 일정하지 않으며 총 이벤트 수에 따라 달라지며 부분적으로는 제어 및 모델링하려는이 효과입니다.

기본 성공률을 평가하기위한 첫 번째 단계. 나는 scipy를 사용하여 그것을 생성함으로써 '관찰 된'데이터의 두 세트를 제공함으로써 상황을 모방하려는 아래의 코드를 준비했다. 그러나 제대로 작동하지 않습니다.
내가 그것을 찾아 낼 것으로 예상 할 것은 :

• loss_lambda_factor 대략 0.1
• total_lambda (및 total_lambda_mu가) 대신 약 120

입니다, 모델이 매우 빠르게 수렴하지만, 예상치 못한에 대답.

• total_lambda 및 total_lambda_mu는 각각 5e5 부근의 날카로운 피크이다. 상기 입력 데이터와 일치하지 않는 번호로 빠르게 수렴, 날카로운 피크 -
• loss_lambda_factor은 (I 인해 10 이하 평판에 게시 할 수 있음) traceplot 상당히 재미가 대략 0

이다. 내가 취하고있는 접근 방식에 근본적으로 잘못된 것이 있는지 궁금합니다. 올바른/예상 된 결과를 얻기 위해 다음 코드를 어떻게 수정해야합니까?

from pymc import Model, Uniform, Normal, Poisson, Metropolis, traceplot 
from pymc import sample 
import scipy.stats 

totalRates = scipy.stats.norm(loc=120, scale=20).rvs(size=10000) 
totalCounts = scipy.stats.poisson.rvs(mu=totalRates) 
successRate = 0.1*totalRates 
successCounts = scipy.stats.poisson.rvs(mu=successRate) 

with Model() as success_model: 
    total_lambda_tau= Uniform('total_lambda_tau', lower=0, upper=100000) 
    total_lambda_mu = Uniform('total_lambda_mu', lower=0, upper=1000000) 
    total_lambda = Normal('total_lambda', mu=total_lambda_mu, tau=total_lambda_tau) 
    total = Poisson('total', mu=total_lambda, observed=totalCounts) 

    loss_lambda_factor = Uniform('loss_lambda_factor', lower=0, upper=1) 
    success_rate = Poisson('success_rate', mu=total_lambda*loss_lambda_factor, observed=successCounts) 

with success_model: 
    step = Metropolis() 
    success_samples = sample(20000, step) #, start) 


plt.figure(figsize=(10, 10)) 
_ = traceplot(success_samples) 

Answer 1

Bayesian MCMC 분석의 함정 (1) 비 집중, (2) 사전, (3) 모델을 제외하고 근본적으로 잘못된 접근 방법은 없습니다.

비 융합 :

이것은 좋은 일이 아니며, 더 명확하게 볼 이유는, 내가 보여주기 위해 traceplot 코드를 변경할 것입니다 :이처럼 보이는 traceplot을 찾을 추적의 두 번째 반, traceplot(success_samples[10000:]) :

이전 : 수렴을위한 주요 도전 과제 중 하나는 베이지안 모델링의 대표적인 함정 인 total_lambda_tau입니다. 이전에 Uniform('total_lambda_tau', lower=0, upper=100000)을 사용하는 것은 매우 유익한 것으로 보일 지 모르지만, 이것은 total_lambda_tau이 큽니다. 예를 들어, 10보다 작은 확률은 .0001입니다. 더 유망하다는 traceplot에서 이전

에

total_lambda_tau= Uniform('total_lambda_tau', lower=0, upper=100) 
total_lambda_mu = Uniform('total_lambda_mu', lower=0, upper=1000) 

결과를 변경 :

이 무엇 나는 traceplot에서 찾아하지만, 더 만족스러운 무언가를 얻기 위해 아직 아니다, "순차적 메트로폴리스 스캔"단계를 사용하는 것이 좋습니다 (PyMC2가 유사한 모델의 기본값 임). 이 허락 된 것 같다 traceplot을 생산

step = pm.CompoundStep([pm.Metropolis([total_lambda_mu]), 
         pm.Metropolis([total_lambda_tau]), 
         pm.Metropolis([total_lambda]), 
         pm.Metropolis([loss_lambda_factor]), 
         ]) 

: 다음과 같이이 작업을 지정할 수 있습니다

이 모델 : @KaiLondenberg는 반응으로 접근하면 total_lambda_tau에 전과로 촬영 한 total_lambda_mu은 표준 방식이 아닙니다. 당신은 광범위하게 변화하는 사건 합계 (1,000 시간과 1,000,000 다음)를 설명하지만, 당신의 모델은 그것이 정상적으로 분배 될 것을 가정합니다.

import pymc as pm, theano.tensor as T 
with Model() as success_model: 
    loss_lambda_rate = pm.Flat('loss_lambda_rate') 
    error = Poisson('error', mu=totalCounts*T.exp(loss_lambda_rate), 
      observed=successCounts) 

나뿐만 아니라 다른 연구 지역 사회에 더 익숙한 것 같다됩니다 다른 방법이 있다는 것을 확신 : 공간 역학에서, 나는 유사한 데이터를 보았다 접근 방식은 더이 같은 모델입니다.

여기에 a notebook collecting up these comments입니다.

Answer 2

모델에 몇 가지 잠재적 인 문제점이 있습니다.

1) 1. 성공 카운트 (오류라고 함)는 포아송이 아닌 이진 (n = 총, p = 손실 _ 감마 계수) 분포를 따라야한다고 생각합니다.

2.) 체인의 시작점은 어디입니까? 순수 Gibbs 샘플링을 사용하지 않는 한 MAP 또는 MLE 구성에서 시작하는 것이 좋습니다. 그렇지 않으면 체인이 번인 (burn-in)에 긴 시간이 걸릴 수 있습니다.

3) total_lambda에 대한 우선적 인 계층 적 선택 (즉, 해당 매개 변수에 대한 두 개의 균일 전치사가있는 정상)은 시작 지점을 현명하게 선택하지 않으면 체인이 수렴하는 데 오랜 시간이 걸릴 수 있습니다 (지점 2와 동일).). 근본적으로 MCMC 체인이 손실 될 수있는 많은 불필요한 자유도를 소개합니다. total_lambda가 비논리적이어야한다는 것을 감안할 때, 적절한 범위 (예 : 0에서 관찰 된 최대 값까지)의 total_lambda보다 먼저 Uniform을 선택합니다.

4. Metropolis Sampler를 사용합니다. 20000 개의 샘플로 충분하지 않을 수 있습니다. 60000을 시도하고 번인으로 첫 번째 20000을 버립니다. Metropolis Sampler는 단계 크기를 조정하는 데 시간이 걸릴 수 있으므로 처음 20000 개 샘플을 사용하여 제안 및 조정을 거부하는 것이 좋습니다. NUTS와 같은 다른 샘플러를 사용해보십시오.