問題タブ [mcmc]

R2jags と CODA を使用して、MCMC チェーンの診断統計を作成していますが、問題があります。次のようにMCMCを実行したい：

エラーは次のとおりです。

RStudio 0.97.551 と R バージョン 3.0.0 (2013-04-03) を使用しています。

どんな助けにも感謝します！

これは私のRスクリプトです：

Xの事前確率(例: X ~ Gaussian) と前方演算子y = f(x)が与えられたとします。さらに、ある実験によってyを観測し、この実験を無期限に繰り返すことができるとします。出力Yは、ガウス (Y ~ ガウス) またはノイズのない (Y ~ デルタ (観察)) と見なされます。

観察されたXに関する主観的な知識の程度を一貫して更新するにはどうすればよいでしょうか? PyMC で次のモデルを試しましたが、何か足りないようです。

事後はxtrueに収束していません。

最初に、私は科学的な数学や統計学の経験が不足していると言わせてください。これは非常によく知られた問題かもしれませんが、どこから始めればよいかわかりません。

f(x1, x2, ..., xn)x'ses を推測し、 の最大値を見つける必要がある関数がありますf。この関数には次のプロパティがあります。

• 合計数またはパラメーターは通常 40 から 60 程度であるため、力ずくのアプローチは不可能です。

• 各 x の可能な値は 0.01 から 2.99 までの範囲です

• 関数は安定しています。つまり、f 値が高いほど、パラメーターの推定が優れていることを意味し、その逆も同様です。

これまでのところ、かなり基本的なメソッドを Python で実装しました。最初にすべてのパラメーターを 1 に設定し、新しい値をランダムに推測して、f が以前よりも高いかどうかを確認します。そうでない場合は、以前の値にロールバックします。10,000 回の反復を伴うループでは、これは何とか機能しているように見えますが、結果はおそらく完璧にはほど遠いものです。

最適なパラメーターの検索を改善する方法についての提案をいただければ幸いです。この問題をグーグルで検索すると、linke MCMC が出てきましたが、それは非常に高度な方法のようで、その方法を理解するのにも多くの時間が必要です。基本的なヒントや概念は、複雑な方法やアルゴリズムよりも役立ちます。

R mcmc コードの C バージョンを作成しています。エラーをスローしているコードの一部は次のとおりです。

blaio.h の関数は、Web サイトの makefile を使用して既にビルドされています。

以下を使用してコンパイルしています。

次のエラーがスローされます。

C にネストされた関数がありません。なぜこのようなエラーが発生するのですか? それがリンクの問題である場合、関数が見つからないようなエラーをスローするべきではありませんか?

編集：解決策を投稿するように求められました。私の最初のコードでは、上記のコードをコンパイルしていませんでした (ただし、別のディレクトリにある古いバージョン)。（うん、それは可能な限り愚かなエラーです）。コードの作業バージョンを再投稿しています。mcmc サンプリング ステップがどのようにコード化され、R からどのように呼び出されるかの例を見る以外に、あまり価値があるとは思いません。ただし、このコードは、グループ化されたベイジアン ラッソに関しては不完全です。コードもあります (blas で使用されるマトリックスを印刷するためのソースが含まれています)。C でのコーディングを開始してちょうど 1 週間ほどです。よりきれいなバージョン (印刷用に別のファイルを配置したり、コマンドをコンパイルするなど) を念頭に置いている場合は、先に進んで変更を加えてください。

この演習で 3 つのことを学びました。

• R から C を呼び出す

• blas を使用して

• gsl を使用して

コードは、後で 4 番目のポイントに対処するように進化します。 - C から R を呼び出す

このモデルでは、一般化された逆ガウス分布からサンプリングを行う必要があります (このコードでは除外されています)。ただし、C でそれを行う簡単な関数は見当たりません。その場合、R の関数を使用する必要があるかもしれません。だから変になる。R はこの C 関数を呼び出し、この C 関数内で他のいくつかの R 関数が呼び出されます。別のアプローチがあれば、コメントをお待ちしております。これが方法である場合は、どうすればよいですか - RInside? もっと簡単ですか？

時間は私にとって最も重要です。少なくとも 10,000 回の mcmc 反復で、約 100 万のデータ ポイントでこのモデルを実行する必要があります。これを5分以内、できれば2分以内に終わらせたいです。

この課題の第 3 段階は、z のサンプリング ステップを並列化することです。

R でランダム データをシミュレートし、C 関数を呼び出す:

技術的な問題があり、自分では解決できないようです。mcmcglmmパッケージで見積もりを実行しました。results$Sol によって、推定事後分布にアクセスできます。class()を適用すると、オブジェクトのクラスが「mcmc」であることがわかります。as.data.frame()を使用すると、他のデータ フレーム (他の多くのデータ フレームを含む 1 つのデータ フレーム) を含むネストされたデータ フレームが生成されます。メイン データ フレーム内のすべてのデータ フレームをrbind()して、すべての事後分布のすべての値と (セカンダリ) データ フレームの名前を行名として持つ 1 つのデータ フレーム (またはベクトル) を生成したいと思います。何か案は？すべてのヒントに感謝します！

更新: Stackoverflow の目的で有用なデータ セットを生成することができませんでした。これらすべてのサンプリング チェーンでは、これらのデータ セットは常に大きすぎます。私を助けたい場合は、次の（例示的な）モデルを実行することを検討してください

result$Sol (mcmc オブジェクト) は、すべてのチェーンが格納される場所です。すべての事後分布のすべての値と変数名を行名として (または追加の文字ベクトルとして重複する行名が許可されていないため) ベクトルにするために、すべてのチェーンを rbind したいと考えています。

私は PyMC に比較的慣れていないので、MCMC サンプラーからの出力に関して簡単な質問があります。MCMC サンプラーによって検出された変数の最も可能性の高い値 (事後変数の最大値) を見つけてほしいと思います。これを行う簡単な方法はありますか？おそらく、MCMC サンプラーによって検出された事後変数の最大値は、PyMC の MAP メソッドによって検出された値とは大幅に異なる可能性があります。

PyMC を提供してくれた開発者に感謝します。それは私の仕事に非常に役立ちます。この質問は、PyMC google グループにも投稿されています。

Rパッケージを使用して現在Rで実行しているJAGS分析のGelmanとRubinの収束診断を計算しようとしていますrjags。

たとえば、パラメーター の収束診断を評価したいと考えていますbeta。これを行うために、ライブラリ coda とコマンドを使用しています。

複数のチェーンを持つout_2MCMCMCMC リスト オブジェクトであるため、エラー メッセージなどのない正しい出力が得られます。ただし、バーンインとして大量の反復を使用したため、反復のサブセットのみ (バーンイン後の部分のみ!) の収束診断を計算したいと思います。

これを行うために、私は試しました：

これにより、次のエラーが発生しました。

したがって、私は試しました：

しかし、驚くべきことに、これにより次のエラーが発生しました。

これは、JAGS 分析から取得したものと同じ MCMC リスト オブジェクトであり、すべての反復の収束診断を評価するときに使用しているものと同じであるため (これは完璧に機能します)、ここに問題はありません。

関数自体は、収束診断の計算でシリーズ (反復) の後半を使用するオプションのみを提供します。私のバーンインフェーズはそれよりも長いので、残念ながらこれでは十分ではありません。

私が見逃しているのは非常に明白なことだと思います。誰か提案やヒントはありますか？

コードが多いため、完全な 2MCMC-JAGS 分析を実行するための R コードは提供しませんでした。上記のコードが問題を十分に示していることを願っています。誰かが以前に同じ問題に遭遇したか、私の構文の間違いを認識した可能性があります。ただし、私の問題を理解するために完全なコードが必要だと思われる場合は、2MCM JAGS 分析を実行するサンプル コードを提供できます。