CでRの拡張したら速すぎて(40〜50倍)吹いたwww - yasuhisa's blog
昨日Gibbs Sampler Algorithmをやってみたわけだが、Rの中でfor文を書いていて必要となるサンプル数が多くなると非常につらくなってくることは目に見えている。しかも、MCMCでは初期値依存となる期間のサンプルを捨てないといけない。そういうわけでじゃんじゃんサンプルを作っても大丈夫なような速度が必要。 Rで速度を上げようと思ったらapplyファミリーを使うとかベクトル単位での処理をするetcが常套手段*1。が、今回は本質的にfor文が必要なケースである。 で、困るわけだがRにはC、C++、fortranを使って拡張する機能がある。詳しくはこの辺に載っている。そういうわけでCのポインタもアドレスも理解していないid:syou6162がRが好きすぎたためにCを書いてみたという感じの内容。 #include <R.h> #include <Rinternals.h> SEXP r
t分布とかX^2分布に従う乱数とか - yasuhisa's blog
初めに言っておくとRならt分布とかX^2分布に従う乱数とかすぐにできるんだけど、あえてそれを正規乱数から生成してみるというやつです。 X^2分布wikipedia:カイ二乗分布より引用すると「を、平均で分散の正規分布に従う、k個の独立なランダム変数とすると、統計量はカイ二乗分布に従う。普通はこれをと書く。カイ二乗分布はkという1個の母数をもつ。」とある。独立性だけ仮定して、同一の分布でなくてもいいのね。忘れていた。というか正規化してあるので、標準正規分布に従う乱数の和と書いてあるのと同値ですね。 なので、ここでは標準正規乱数の和を用いてX^2分布に従う乱数を生成してみましょう。以下では標準正規乱数を1000個作っています。 rnorm(1000,0,1) 標準正規乱数の和がX^2分布に従う乱数ということなので、以下の操作で自由度1000のX^2分布に従う乱数が一つ生成できました*1。やった
1
リリース、障害情報などのサービスのお知らせ
最新の人気エントリーの配信
処理を実行中です
j次のブックマーク
k前のブックマーク
lあとで読む
eコメント一覧を開く
oページを開く
http://cran.r-project.org/doc/contrib/R-and-octave.txt
