Rによる確率分布からのサンプリング

　実験コードが不可思議な挙動を示すので、小さなデータセットを作って意図通りの挙動を示すかどうか確認してみることになった。なぜ修論提出後にこんな事をしているのかというと、実装が間に合わなかったからである。これについてはまたいつか。
　LDAに従った文書生成というのも考えたが、そこまでする意味もないだろうということで、もっと単純に、トピック毎に多項分布を用意し、そこからサンプリングした文書を混ぜて使うことにした。
　さて、多項分布からのサンプリングなんてコード、どうやって書けばいいのか。一様乱数で良ければいろんな乱数ジェネレータがすぐに見付かる。正規分布からのサンプリングもまぁ、なんか教科書で見たような気がする。一様乱数があれば正規分布もなんとかなるだろう。しかし、多項分布はどうすればいいのか。
　と、ここで無理に自分でコードを書いていたらそれだけで1ヶ月とかかかりそうだったので、Rを使ってサンプリングをしてみることにした。Rを使うのは初めてである。r-baseとかいうパッケージをインストールしてみた。起動コマンドはRらしい。大文字で。
　rnorm(欲しい乱数の個数,m=平均値,sd=標準偏差)で、正規分布からのサンプリングができるそうだ。第一引数個分の要素を含んだベクトルが返ってくる。変数に束縛したい場合は<-を使ってvalues <- rnorm(10, m=1, sd=1)みたいにする。多項分布からのサンプリングには同じようにrmultinomという関数（？）を使う。rmultinom(個数, 多項分布のパラメータ数, パラメータの確率分布)みたいな感じでサンプリングができる。行列が返ってくるので、後はそれを保存すればよろしい。　
　保存にはwrite.tableという関数を使う。（ピリオドも関数名に含められるようだ。）オプションで行列名を保存するとかしないとかが選べる。（col.name=Fとかrow.name=Fとかがそれ。）名前つきパラメータって、使う方にしてみれば便利だなぁ。
行列を転置する場合の関数はt。ベクトルの連結にはc。素晴らしく短い。Schemeならmatrix-transposeとかvector-concatenateとか、そんな名前を付けてるところなのではないかと思う。
　こうして、以下の3行のコードで、多項分布からのサンプリングと、そのファイルへの保存が出来た。

topicA <- c(rnorm(10,m=0.5,sd=0.1), rnorm(10,m=1.0,sd=0.1), rnorm(10,m=0.2,sd=0.1))
topicA_data <- rmultinom(225, size = 200, topicA)
write.table(t(topicA_data),"topicA_data",quote=F,col.name=F,row.name=F,append=F)

topicAを正規化しておく必要は無い。ちょっとしたことだが、気が利いてて嬉しい。
これだけの事を調べるのに、半日ほど潰した。検索しにくい名前だったのでなかなか辛かった。
試してないが、MCMCpackというのを使うと、Dirichlet分布からのサンプリングも行えるそうだ。