問題タブ [mcmc]

2 つのパラメーターを推測するために、PyMC3 で動的システムのモデルをまとめようとしています。モデルは、疫学で一般的に使用される基本的な SIR です。

dS/dt = - r0 * g * S * I

dI/dt = g * I ( r * S - 1 )

ここで、r0 と g は推論されるパラメーターです。これまでのところ、私はまったく遠くに行くことができません。このようなマルコフ連鎖をまとめた唯一の例では、再帰が深すぎるというエラーが発生します。これが私のコード例です。

どんな助けでも大歓迎です。ありがとう ！

pymc2 を pymc3 に使用したいくつかの計算を更新しています。モデルに離散確率変数がある場合、サンプラーの動作に問題があります。例として、pymc2 を使用した次のモデルを考えてみましょう。

これは何かを代表するものではなく、観測されていない変数の 1 つが離散しているモデルにすぎません。このモデルを pymc2 でサンプリングすると、次の結果が得られます。

しかし、PYMC3 で同じことを試すと、次のようになります。

変数 A がまったくサンプリングされていないようです。pymc3で使用されているサンプリング方法についてはあまり読んでいませんでしたが、特に連続モデルを対象としているように見えました。これは、モデルで観測されていない離散変数を除外することを意味しますか、それとも私がやろうとしていることを行う方法はありますか?

(X,Y) ポイントのランダムなコレクションがあるとします。

そして、単純な線形回帰に適合すること：

y_est[0]、y_est[1]、 ..の分布または統計を取得するにはどうすればよいですか (これらの変数は、各入力値の推定値y_est[2]に対応することに注意してください。yx

データから単純な線形関数とガンマ分布ノイズ項のパラメーターを推定したいと思います。(注: これはhttps://stats.stackexchange.com/questions/88676/regression-with-unidirectional-noiseのフォローアップの質問ですが、単純化され、より実装固有です)。次のように生成された観測データがあるとします。

次のようになります。

次のように、pymc を使用してこれらのパラメーターを推定してみました。

ただし、これにより、真の値からかけ離れた見積もりが得られます。

• 切片: 真: 1.000、推定: 3.281
• 勾配: 真: 2.000、推定: -3.400

私は何か間違ったことをしていますか？

output = RBugsfit(..., coda=T, ...)4 つのパラメーターの事後分布のサンプルとそれらのサンプル事後平均を含む mcmc.list オブジェクトを出力します。を使用しsummary()てサンプルの事後平均を確認できますが、サンプルの事後平均をoutputプログラムの変数に取得する方法を知りたいですか? ありがとう！

ステップごとにサンプリングするときにサンプラーが非常に遅い理由を知りたいです。たとえば、次のように実行すると:

サンプリングは高速です。ただし、実行すると：

サンプリングは遅くなります (サンプリングが多いほど遅くなります)。

なぜ私がこれを尋ねているのか疑問に思っているなら..まあ、サンプラーの各ステップの後に変数の値に対していくつかの操作を実行したいので、段階的なサンプリングが必要です。

それをスピードアップする方法はありますか？

前もって感謝します！

- - - - - - - - - 編集 - - - - - - - - - - - - - - - - ------------------------------

ここで、特定の問題をより詳細に提示します。

競合する 2 つのモデルがあり、それらは 2 つの間の「スイッチ」として機能するカテゴリ変数を持つより大きなモデルの一部です。

このおもちゃの例では、観測されたベクトル 'Y' があり、これはポアソン分布または幾何分布で説明できます。カテゴリ変数 'switch_model' は、= 0 の場合は幾何モデルを選択し、=1 の場合はポアソン モデルを選択します。

各サンプルの後、switch_model が幾何学的モデルを選択した場合、ポアソン モデルの変数は更新されないようにします。これは、可能性に影響を与えていないため、単に漂流しているだけだからです。switch_model がポアソン モデルを選択する場合は、その逆です。

基本的に、各ステップで行うことは、選択されていないモデルの値を手動で 1 ステップ戻すことによって「変更」することです。

私の説明とコメント付きのコードが十分に明確になることを願っています。詳細が必要な場合はお知らせください。

モデルのすべてのパラメーターを推定したい多項モデル (4 つのモダリティがあります) を使用しています (3 つの最初のモダリティのそれぞれに 2 つのパラメーターがあるため、6 つのパラメーター)。私はそれを行うことができました。私の問題は、パラメーターに制約を追加して、最初のモダリティの確率が共変量で減少するようにしたいということxです。

したがって、次のものが必要であると計算します。

runjagsJAGSを使ったRのパッケージを使っています。各サンプルについて、この制約が尊重されることを望みます。それは可能ですか？

通常は C/C++ ライブラリを使用し、各反復で条件をチェックし、条件が満たされない場合はこの反復のすべてのパラメーターを再サンプリングします。しかし、私は JAGS でそれを行う方法を見つけることができません (また、コラボレーションの理由から、C を使用できません)。

現在のコードは次のとおりです。

よろしくお願いいたします。

一番

Haskell で Metropolis アルゴリズムに基づいて MCMC プログラムを作成しようとしていますが、確率分布からのサンプリング (疑似乱数の生成) とプログラムの構造化に問題があります。今のところ、IO を処理する複雑さに対処するよりも、ハードコードされたシードを持つジェネレーターを使用することに満足しています。

It seems that I should use the state monad to keep track of the random generator state, the previous markov chain state, chi squared value and the acceptance count between each step of the algorithm and then finally collect all the markov chain states and the last acceptance count. Is this the best/idiomatic way to do this? and if so what should the layout of the program should be (i.e. type signatures of the proposal function and the metropolis step function etc.).

I have seen some sample programs that deal with random numbers where a list of random numbers of a particular length is generated from some kind of a probability monad and then threaded through some simple functions to perform the computation. I really would like to avoid this inside out form of the program if it possible.

Edit: Temporarily removed WIP code.