10

パラメータを推定したい観測データがあり、PYMC3を試す良い機会だと思いました。

私のデータは一連のレコードとして構造化されています。各レコードには、固定された 1 時間の期間に関連する観測のペアが含まれています。1 つの観測値は、特定の時間内に発生したイベントの総数です。もう 1 つの観察結果は、その期間内の成功数です。たとえば、データ ポイントは、特定の 1 時間の間に合計 1000 件のイベントがあり、そのうちの 1000 件のうち 100 件が成功したことを示す場合があります。別の期間では、合計で 1000000 のイベントが発生し、そのうち 120000 が成功したとします。観測値の分散は一定ではなく、イベントの総数に依存します。私が制御およびモデル化したいのは、この効果の一部です。

これを行うための最初のステップは、基本的な成功率を推定することです。以下のコードは、scipy を使用して生成することで 2 セットの「観察された」データを提供することにより、状況を模倣することを目的としています。ただし、正しく動作しません。
私が期待しているのは次のとおりです。

  • loss_lambda_factor はおよそ 0.1 です
  • total_lambda (および total_lambda_mu) はおよそ 120 です。

代わりに、モデルは非常に迅速に収束しますが、予想外の答えになります。

  • total_lambda と total_lambda_mu は、それぞれ 5e5 付近の鋭いピークです。
  • loss_lambda_factor はほぼ 0 です。

traceplot (レピュテ​​ーションが 10 未満であるため投稿できません) はかなり面白くありません。収束が速く、入力データに対応しない数値で鋭いピークが発生します。私が取っているアプローチに根本的な問題があるかどうか、私は興味があります. 正しい/期待される結果を得るには、次のコードをどのように変更する必要がありますか?

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)
4

2 に答える 2

33

ベイジアン MCMC 分析の落とし穴 (1) 非収束、(2) 事前確率、(3) モデルを除いて、アプローチに根本的な問題はありません。

非収束: 次のような traceplot を見つけました。

バーニンを含む traceplot

これは良いことではありません。その理由をより明確に理解するために、トレースプロット コードを変更して、トレースの後半のみを表示しますtraceplot(success_samples[10000:])。

バーニンを削除した traceplot

事前確率 : 収束の主要な課題の 1 つは、事前確率 ontotal_lambda_tauです。これは、ベイジアン モデリングの典型的な落とし穴です。before を使用することはまったく有益ではないように見えるかもしれませんが、これの効果は、が大きいUniform('total_lambda_tau', lower=0, upper=100000)ことを確信していると言うことです。total_lambda_tauたとえば、10 未満である確率は .0001 です。以前の変更

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

より有望な traceplot が得られます。

事前分布が異なる traceplot

ただし、これはまだトレースプロットで探しているものではありません。より満足のいくものを得るには、「順次スキャン メトロポリス」ステップを使用することをお勧めします (これは、類似モデルに対して PyMC2 がデフォルトで使用するものです)。これは次のように指定できます。

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

これにより、許容できると思われるトレースプロットが生成されます。

シーケンシャル スキャン メトロポリスを使用した traceplot

モデル: @KaiLondenberg が回答したように、事前に取ったアプローチは標準的なアプローチではありませんtotal_lambda_tau。total_lambda_muさまざまなイベントの合計 (1 時間で 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)

他の研究コミュニティでも、より馴染みのある方法が他にもあると確信しています。

これは、これらのコメントをまとめたノートです。

于 2014-06-17T19:21:30.543 に答える
1

モデルには潜在的な問題がいくつかあります。

1.) 成功数 (エラーと呼ばれますか?) は、ポアソン分布ではなく、2 項 (n=total,p=loss_lambda_factor) 分布に従うべきだと思います。

2.) チェーンはどこから始まりますか? 純粋なギブス サンプリングを使用しない限り、MAP または MLE 構成から開始することは理にかなっています。そうしないと、チェーンがバーンインするのに長い時間がかかる可能性があり、それがここで起こっている可能性があります。

3.) total_lambda の階層的な事前確率 (つまり、これらのパラメーターに 2 つの均一な事前確率を持つ通常) を選択すると、(ポイント 2. のように) 開始点を賢く選択しない限り、チェーンが収束するのに長い時間がかかることが保証されます。基本的に、MCMC チェーンが迷子になるために、多くの不必要な自由度が導入されます。

4.) メトロポリス サンプラーを使用します。20000 サンプルでは十分ではない可能性があります。60000 を試し、最初の 20000 をバーンインとして破棄します。Metropolis Sampler はステップ サイズの調整に時間がかかることがあるため、最初の 20000 サンプルを主に提案の拒否と調整に費やした可能性があります。NUTS などの他のサンプラーを試してください。

于 2014-06-17T08:07:27.380 に答える