Hatena::ブログ(Diary)

My Life as a Mock Quant このページをアンテナに追加 RSSフィード Twitter

2014-10-11

逆ガンマ分布(事後分布)〜正規分布(尤度関数)×逆ガンマ分布(事前分布)

| 10:34 | 逆ガンマ分布(事後分布)〜正規分布(尤度関数)×逆ガンマ分布(事前分布)を含むブックマーク

一般に、ある確率測度¥mathbb{P}の下で、逆ガンマ分布に従う確率変数X確率密度関数 f(x; ¥alpha, ¥beta)

 f(x; ¥alpha, ¥beta) = ¥frac{¥beta^¥alpha}{¥Gamma(¥alpha)} x^{-¥alpha - 1}¥exp¥left¥{ -¥frac{¥beta}{x}¥right¥}

である。

この時、確率変数X分散に持つ正規分布N(¥mu, X)を尤度関数とし、逆ガンマ分布 f(x; ¥alpha_0, ¥beta_0)を事前分布とすると、

確率変数X事後確率密度関数 ¥tilde{f}(x; ¥alpha, ¥beta)

 ¥tilde{f}(x; ¥alpha, ¥beta) ¥propto ¥frac{1}{¥sqrt{2¥pi}x} ¥exp¥left¥{ -¥frac{(y-¥mu)^2}{2x}¥right¥}  ¥frac{¥beta_0^{¥alpha_0}}{¥Gamma(¥alpha_0)} x^{-¥alpha_0 - 1}¥exp¥left¥{-¥frac{¥beta_0}{x}¥right¥} ¥propto ¥frac{1}{x} ¥exp¥left¥{ -¥frac{(y-¥mu)^2}{2x}¥right¥}  x^{-¥alpha_0 - 1}¥exp¥left¥{-¥frac{¥beta_0}{x}¥right¥} ¥propto x^{-¥alpha_0 - 2} ¥exp¥left¥{ -¥frac{1}{x}¥left( ¥beta_0 + ¥frac{1}{2}(y-¥mu)^2¥right)¥right¥}

従って、

 ¥alpha = ¥alpha_0 +1

 ¥beta = ¥beta_0 + ¥frac{1}{2}(y-¥mu)^2

となり、 ¥tilde{f}(x; ¥alpha, ¥beta)パラメーターが決定され、かつ、その確率分布自体は逆ガンマ分布となる。すなわち自然共役である。

参考

2014-10-10

Rにおける日付(Date型)のデフォルトの原点は1970年1月1日だぞっと

| 21:09 | Rにおける日付(Date型)のデフォルトの原点は1970年1月1日だぞっとを含むブックマーク

前振り

数値を日付に変換する時には基準となる日付が必要なわけですが、Rの場合、デフォルトだと「1970年1月1日」になってるようだ。

この話のモチベーションとしては別に俺的には日付のまま処理したいんだけど、lapplyの結果を行列にしたりデータフレームにしたりしてるといつの間にか日付が数値になってて、どうやって直そうかとアワアワするのをなくすことだ。

本題

例えば今日の日付をGETしておく

> today <- Sys.Date()
> today
[1] "2014-10-10"

そして、それを数値にしてみる。

> today_as_number <- as.numeric(today)
> today_as_number
[1] 16353

これを日付に戻そうとas.Date関数を噛ませると、originを指定してないと怒られる。

> as.Date(today_as_number)
Error in as.Date.numeric(today_as_number) : 
   'origin' を指定しなければなりません 

ここで、originとして1970年1月1日を指定すると元の日付に戻るので、これで正しそうだぞとわかる。

> as.Date(today_as_number, origin="1970-01-01")
[1] "2014-10-10"

lubridate使えばこんなことにはならないのかなぁ。

ここに答えがあったんや

- [Rd] as.Date.numeric origin default

2014-09-28

正規分布(事後分布)〜正規分布(尤度関数)×正規分布(事前分布)

| 08:23 | 正規分布(事後分布)〜正規分布(尤度関数)×正規分布(事前分布)を含むブックマーク

一般に、ある確率測度¥mathbb{P}の下で、正規分布に従う確率変数Xに対し、

  • 尤度関数X ¥sim ¥mathcal{N}(¥mu_1, ¥sigma_1)
  • 事前分布:X ¥sim ¥mathcal{N}(¥mu_0, ¥sigma_0)

であることを仮定したとき、事前分布が共役事前分布となっていることから、事後分布も正規分布 ¥mathcal{N}(¥hat{¥mu}, ¥hat{¥sigma})となり、

  •  ¥frac{1}{¥hat{¥sigma}^2} = ¥frac{1}{¥sigma_1^2} + ¥frac{1}{¥sigma_0^2}
  •  ¥hat{¥mu} = ¥frac{¥hat{¥sigma}^2}{¥sigma_1^2} ¥mu_1 + ¥frac{¥hat{¥sigma}^2}{¥sigma_0^2} ¥mu_0

が成立する。ただしここで、¥mathcal{N}(¥mu, ¥sigma)は平均¥mu, 分散¥sigmaを持つ正規分布を表す。従って、データが N ¥in ¥mathbb{Z}個存在する場合には、上式を再帰的に用いることで

  •  ¥frac{1}{¥hat{¥sigma}^2} = ¥sum_{i=1}^{N} ¥frac{1}{¥sigma_i^2} + ¥frac{1}{¥sigma_0^2}
  •  ¥hat{¥mu} = ¥sum_{i=1}^{N} ¥frac{¥hat{¥sigma}^2}{¥sigma_i^2} ¥mu_i + ¥frac{¥hat{¥sigma}^2}{¥sigma_0^2} ¥mu_0

が成立する。

2014-09-23

foreachパッケージで、並列化する/しないを綺麗に書き分ける

| 08:14 | foreachパッケージで、並列化する/しないを綺麗に書き分けるを含むブックマーク

で、「これからRで並列化するならforeachパッケージ使うだろ」という話を書いたが、その続き的な?もの。

並列化する/しないを、どう美しく書き分けるにはどうしたらよいかというお話。

たとえば、自作の関数かなんかで、並列化する/しないを書き分けようとすると、引数かなんかにparallel引数を用意しておいて、以下のように書いてしまうんじゃかろうか。

parallel=TRUE
if(parallel){
  foreach(i=1:3) %dopar% {
    rnorm(i)
  }
}else{
  foreach(i=1:3) %do% {
    rnorm(i)
  }  
}

これはdo/dopar部分だけが違うだけであって、冗長なコードだ。Rだと、最近流行りのパイプ一族系ポスト

を見てもわかるように、%で挟まれた奴らは実はこれ関数なのです。なので、以下のように書くことで、上のコードは

parallel=TRUE
doop <- ifelse(parallel, `%do%`, `%dopar%`)
doop(foreach(i=1:3), {rnorm(i)})

とまで簡略化される。こっちのほうがスマートに(私には)見えるので、これでいく。ここのコードではワーカーの立ち上げ・停止は書いていないことに注意。立ち上げたワーカー(スレッド、というかRScript.exeプロセス)を正しく停止させるためには以下を見ておくとよい。

これを受けて、pipeR三銃士

である俺はpipeRパッケージを活用して以下のように書いておくことにした。

library(pipeR)
parallel <- TRUE
doop <- if(parallel){detectCores() %>>% (cl = makeCluster(.)) %>>% registerDoParallel; `%dopar%`}else{`%do%`}
doop(foreach(i=1:3), {rnorm(i)})
if(parallel){stopCluster(cl)}

さようなら可読性。そしてさらにコメント欄で貰ったアドバイスをもとにすると

library(pipeR)
parallel <- TRUE
"%doop%" <- if(parallel){detectCores() %>>% (cl = makeCluster(.)) %>>% registerDoParallel; `%dopar%`}else{`%do%`}
foreach(i=1:3) %doop% {
  rnorm(i)
}
if(parallel){stopCluster(cl)}

とできる。これがスマートかなぁ。

通りすがり通りすがり 2014/09/23 22:14 いっそdoopも演算子にしてしまえば良いのではないでしょうか。

"%doop%" <- ifelse(parallel, `%do%`, `%dopar%`)

teramonagiteramonagi 2014/09/24 06:17 それだ!

2014-09-20

foreachパッケージで並列化する時、現在の"環境"にない変数・関数は、明示的に.export引数にて指定しなければならない

| 15:39 | foreachパッケージで並列化する時、現在の"環境"にない変数・関数は、明示的に.export引数にて指定しなければならないを含むブックマーク

foreachパッケージについて

遅い遅いと巷で噂のR君も並列化すればそれなりにパフォーマンスが期待できるわけです。Rにおいてどうやって並列化をするのかというと、foreachパッケージ&関数で並列化するのが手っ取り早くて、これはざっくりでいうと

foreach(i in 1:3) %dopar%{
    #ここに処理
}

と書くだけで処理の並列化をすることが可能です。ここで%dopar%と書いている箇所を%do%に直せば並列化させないとして処理されるのも良いところだ。この裏側の仕組みとしてはforeachパッケージが各種並列化パッケージであるsnow/multicore/parallel/Rmpiへのフロントエンドとして機能していて、裏側ではそれらのパッケージに処理をブン投げる形になっているそうだ。詳しくは

を読むとよい。上述の書籍を読むに、このforeachパッケージとparallelパッケージを組み合わせた並列化がこれからのデファクトスタンダードになりそうなんで、ここではparallel+foreachパッケージでやってみる。というわけでパッケージのロードをしておこう。

library(foreach)
library(doParallel)

さらに並列化の準備として

registerDoParallel(detectCores())

も実行しておく。これはdetectCores関数で現在のマシンのコア数を取得し、それらを並列環境としてregisterDoParallel関数で登録しますよということだ。

俺がやりたかったことと、ハマったこと

適当な乱数生成をする関数を以下のように2つ用意する。f1のほうが並列化した版で、f2の方が並列化していない版に相当する。

f1 <- function()
{
  foreach(i=1:3) %dopar% {rnorm(i)}
}
f2 <- function()
{
  foreach(i=1:3) %do% {g(i)}
}

これを実行すると、正しく求めたい結果が返ってきて

> f1()
[[1]]
[1] 0.4182321

[[2]]
[1] 0.8150335 1.2575948

[[3]]
[1]  0.5326094 -1.3211135  1.7029272

> f2()
[[1]]
[1] -0.7583267

[[2]]
[1]  0.5325022 -0.4120185

[[3]]
[1]  1.58705166 -0.58542707  0.03252658

となる。一方、実際に並列化させたい処理なんかは、自作の関数を用いることも多々あるわけで、ここではその自作関数の例として関数gを定義して、さきほどと同じことをやろうとする・・・と・・・

g <- function(i){rnorm(i)}
f3 <- function()
{
  foreach(i=1:3) %dopar% {g(i)}
}

↓実行すると・・・

> f3()
Error in { : 
  task 1 failed - " 関数 "g" を見つけることができませんでした "
Called from: foreach(i = 1:3) %dopar% {
    g(i)
}

となって怒られる。これはどうも処理を並列化する際に、並列化するスレッドに対して適切な環境*1が渡っていないからのように見える。

解決策

この問題を解決するには明示的にgを並列計算する側の環境上にもexport(輸出)しますよーというのを、その名の通りの.export引数で明示すればいい。

f4 <- function()
{
  foreach(i=1:3, .export="g") %dopar% {g(i)}
}

こうすることで、ちゃんと結果を返すようになる。

> f4()
[[1]]
[1] -0.07732738

[[2]]
[1] -0.6724723  0.5245355

[[3]]
[1]  0.2095755 -0.2227492 -0.4902599

として実行することができる。もうちょっといろいろいじってみたところ、どうもforeachパッケージを用いると

  • foreachに実際に処理を移譲するコードを描いた環境の情報(ここだと関数f4内部で定義されている変数関数のみ)

しか並列計算させる環境の上には持っていってくれないということだ。同様のことはパッケージについても言えて、自作の並列化処理において、例えばXXX, YYYパッケージを用いているならば、foreach関数引数として

foreach(i=1:10, .packages=c('XXX', 'YYY'))

というようにそれらを「明示的に使いますよ」という記述が必要になる。

もうちょっと楽に

今は関数g一個を並列化させる環境上にexportさせればよかったが、もっと多数の変数関数をexportする場合に、いちいちこれを全部書きあげるってのはとても面倒だ。なので、こういう面倒なことをするのが億劫な人、特に俺は

ls(envir=parent.frame())

を.exportの引数に入れておくことにしたい。こうすると親の環境(メインの環境だと思われるところ)の情報(関数とか変数)を全部持って行ってくれるので、楽。

f4 <- function()
{
  foreach(i=1:3, .export=ls(envir=parent.frame())) %dopar% {g(i)}
}

これでもちゃんと動作することは確認できて

> f4()
[[1]]
[1] -0.04673488

[[2]]
[1] -0.3213969 -1.3645549

[[3]]
[1] -1.9164478  0.3883101  2.2534236

となる。実際にls(envir=parent.frame())がどんな出力を返しているのかは上でいう関数f4の中でprint(ls(envir=parent.frame()))というように出力させてみるといい。出力文字列の中にgという文字がみつかるはずだ。




・・・みたいな話も上述の並列化本に書いてあった・・・ちゃんと読んでから並列化すればよかったorz

*1:自分で作ったRの変数関数の集合というか、定義したもろもろのものと考えればOK

774774 2014/09/20 17:06 「上述の書籍」とは?

teramonagiteramonagi 2014/09/21 07:59 「Rによるハイパフォーマンスコンピューティング」です!