Stick breaking process in stan
The BUGS bookの293 pageに書かれていたStick breaking processをstanで実装、コンパイルができましたが、Errorの発生により評価が阻まれています。
82個の銀河の銀河系からの相対速度の分布がテストデータとして使われていますが、入手が難しそうなので手で作成、試行します。
code
実質的クラスタ数(BUGS codeのK)や推定された密度分布(BUGS codeのdens)はR codeで計算することとします。
Error message
Invalid value of pi: Error in function validate transformed params N4stan5agrad3varE: pi
pi is not a valid simplex. The element at 5 is -0.719628:0, but should be greater than or equal to 0
Rejecting proposed initial value with zero density.Initialization failed after 100 attempts. Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
error occurred during calling the sampler; sampling not done
sampling中にpiが負の値になってしまっている場合があるようです。
Reference
データの元論文
http://www.cs.berkeley.edu/~jordan/courses/281B-spring04/readings/escobar-west.pdf
http://www2.imm.dtu.dk/courses/02443/projects/Roeder_JASA_1990.pdf
incremental_log_probについて
http://heartruptcy.blog.fc2.com/blog-entry-141.html
http://heartruptcy.blog.fc2.com/blog-entry-155.html
http://qiita.com/hoxo_m/items/56a63faf29e1af4db27d
http://xiangze.hatenablog.com/entry/2013/12/19/013557
StanのCategorical分布をindexとして使えない部分はモデリング言語としての可読性を著しく下げていますが、計算速度の点からは余計な分布を計算しなくて良いので、故あってのことといえるかもしれません。
@hoxo_m @teramonagi @kos59125 例えばstan 2.4.0のマニュアルのp.55 Analytic Posteriorでは共役事前分布を使って式変形することで高速化が可能になっています。その他の例→ http://t.co/ILeQReazjf
— berobero (@berobero11) August 2, 2014
https://groups.google.com/forum/#!msg/stan-users/LK0kGNdyGdE/ajfjMe4C1ngJ
ノンパラメトリックベイズモデルについて
http://cran.r-project.org/web/packages/DPpackage/index.html
パッケージユーザーのための機械学習(9):混合ディリクレ過程
RにはDPpackageがあるそうです。The BUGS book 293pageでも紹介されています。
PyMC 3でのstick-breaking processの実装だそうです。
Categorical variable always zero · Issue #604 · pymc-devs/pymc · GitHub
CRPの方はstanでの実装は難しそうです。