【2021-06-15 更新】本記事の内容は時代遅れです。 宋財泫先生のチュートリアルを見ましょう。

purrrと関数

みんな大好きtidyverseの一部を成す purrrというライブラリーを使って 「関数のリスト」にリストを投げる方法を紹介します。 まずは「なんでそんな状況が起きうるのか」という話からはじめます。 ちなみに”purr”というのは猫のゴロゴロのことで、リーナスが好む雑音の一つです。

さて、確率分布の一つにポアソン分布というのがあります。 ポアソン分布は平均(lambda)を引数にカウントデータの分布を説明してくれます。 データに含まれる個々の値をとり、 それぞれが自身（平均lambdaで作られたポアソン分布）に属す確率を返す、という具合です。

# ラムダをとってポアソン分布（関数）を返す。lはlambdaのl。
pois <- function(l) function(k) l^k*exp(-l)/factorial(k)

まずはこの関数がどう動くのかを見てみます。 久保先生のサイトからデータを引っ張ってきて 「平均3.56のポアソン分布」にそれぞれのデータが属する可能性を見てみます。 50個のデータのそれぞれが「平均3.56のポアソン分布」に属する確率が出力されてますね。 

load(url("http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/distribution/data.RData"))
pois(3.56)(data)
## [1] 0.18021114 0.18021114 0.19032700 0.08040427 0.19032700 0.13551282
## [7] 0.18021114 0.21385056 0.10124222 0.18021114 0.02843882 0.19032700
##[13] 0.21385056 0.21385056 0.21385056 0.21385056 0.19032700 0.18021114
##[19] 0.04089132 0.18021114 0.19032700 0.21385056 0.21385056 0.21385056
##[25] 0.19032700 0.21385056 0.04089132 0.13551282 0.21385056 0.10124222
##[31] 0.04089132 0.08040427 0.19032700 0.08040427 0.13551282 0.18021114
##[37] 0.19032700 0.04089132 0.18021114 0.18021114 0.08040427 0.18021114
##[43] 0.19032700 0.13551282 0.19032700 0.13551282 0.10124222 0.21385056
##[49] 0.18021114 0.21385056

このような確率分布のリストを3.56だけではなく、色々なパラメターで作りたいケースを考えます。 モチベーションとしては「与えられたデータがどの平均のポアソン分布に属するのか」 の推定があります。 つまり、パラメター(データの平均)となりそうな値の周辺のリストを 平均の候補としてポアソン分布へ適用します。 するとポアソン分布のリストができあがるわけです。

関数のリストを作成する

先程は勝手にパラメターを3.56と決めました。 それ以外のパラメターを持つポアソン分布ではどの程度あてはまるのか、 というのが次に見たいポイントです。

そこでlambdasという変数に2から5の0.1刻みのatomicなlistを作ります。 名前の通り、lambdaが複数(31個)入っています。

lambdas = seq(2, 5, 0.1)

そしてそれぞれの値でポアソン分布を作ります。 そのためにはpurrrパッケージのmap 関数を使います。 このmap関数は、lambdasの中にあるlambdaひとつひとつを 第二引数にあるpois、つまり先ほど作った関数に適用する、というものです。

library(purrr)
map(lambdas,pois)
## [[1]]
## function (k) 
## l^k * exp(-l)/factorial(k)
## <bytecode: 0x36411a8>
## <environment: 0x33f9a78>
## :
## :
## [[31]]
## function (k) 
## l^k * exp(-l)/factorial(k)
## <bytecode: 0x36411a8>
## <environment: 0x3649ae0>

これで31個のポアソン分布、関数が生成されるわけです。 それぞれがdataに対してどれくらいの当てはまりの良さを持っているかは、 それぞれがdataを引数にとってdataの各要素がそれぞれの分布に属す確率をもとめ、 掛けあわせればいいのですね。 しかし実際にすると色々と問題があるので、ログをとって足す、という作業をします。

library(magrittr)
# lambdas = seq(2, 5, 0.1)
lambdas %>% map(pois) %>%
     # map(pois)で作った関数のリストを"."として渡す
     map(., ~map_dbl(data, .)) %>%
     map(log) %>% map_dbl(sum)
##[1] -121.88118 -118.19653 -114.91597 -112.00356 -109.42794 -107.16163
##[7] -105.18034 -103.46256 -101.98912 -100.74287  -99.70839  -98.87180
##[13]  -98.22054  -97.74318  -97.42935  -97.26957  -97.25516  -97.37814
##[19]  -97.63119  -98.00755  -98.50098  -99.10570  -99.81633 -100.62791
##[25] -101.53577 -102.53560 -103.62336 -104.79525 -106.04775 -107.37751
##[31] -108.78143

関数のリストを実行する

以下の部分を詳しく説明すると以下の流れになります

  1. 平均のリストをpoisに投げて(%>%)分布(関数)のリストを作る
  2. ポアソンのリストそれぞれを~map_dbl(data, .)の.に置く
  3. するとdataをポアソンのリストに投げれる
lambdas %>% map(pois) %>%
     # map(pois)で作った関数のリストを"."として渡す
     map(., ~map_dbl(data, .))

それを更にlogに投げてsumに投げるんですね。 ポアソン分布を使って尤度を図ろうとすると 定義と実装がずれてしまって嬉しくない場合があります。 purrrとか使って関数のキャッチボールができると より自然な書き方ができて非常に嬉しいなと思う日々です。

参考になりそうな文献