この企画は、雑誌や教科書にでているグラフをSageで再現し、 グラフの意味を理解すると共にSageの使い方をマスターすることを目的としています。
今回は、データ解析のための統計モデリング入門 (以下、久保本と書きます) の第2章の例題をSageを使って再現してみます。久保氏のブログは、WinBUGSの使い方などでよく拝見していましたが、 「データ解析のための統計モデリング入門」は自然科学を学ぶすべての人に薦めたい一冊です。まさに目から鱗の固まりです。
数式処理システムSageのノートブックは、計算結果を表示するだけではなく、実際に動かすことができるの大きな特徴です。 この機会にSageを分析に活用してみてはいかがでしょう。
最初に必要なライブラリーやパッケージをロードしておきます。新しくなったRUtil.pyも使います。
|
2章の例題に使われているデータは、架空の植物の第i番目の個体の種子数$y_i$を扱っています。 本の図2.1を以下に引用します。
2章のデータは、ネット公開されており、本のサポートサイトからダウンロードできます。 ここでは、DATAディレクトリにdata.RDataをアップし、このノートブックで参照できるようにしました。
データをロードしたら、summary関数とvar関数を使ってデータの素性を確認します。 ここで、確認することは、以下の3項目です。
Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00 2.00 3.00 3.56 4.75 7.00 [1] 2.986122 Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00 2.00 3.00 3.56 4.75 7.00 [1] 2.986122 |
データの分布を見るには、Rのtable関数を使って度数分布を計算すると便利です。またグラフに表示するときには、hist関数を使用します。
data 0 1 2 3 4 5 6 7 1 3 11 12 10 5 4 4 data 0 1 2 3 4 5 6 7 1 3 11 12 10 5 4 4 |
|
Rのライブラリはとても豊富で実績もありますが、個人的には関数の定義が分かりずらいのであまり好きではありません。
そこで、SageにPandasライブラリをインストールして、これを使って同じような処理をしてみました。
Rから配列データを取り出すには、sageobj関数を使います。
|
Rから取り込んだ2章のデータをPandasのデータフレームにします。
|
データフレームにしたorgDataからRのsummaryと同様の情報をdescribeを使って得ることができます。
data count 50.00000 mean 3.56000 std 1.72804 min 0.00000 25% 2.00000 50% 3.00000 75% 4.75000 max 7.00000 data count 50.00000 mean 3.56000 std 1.72804 min 0.00000 25% 2.00000 50% 3.00000 75% 4.75000 max 7.00000 |
Sageにはヒストグラムをプロットする関数がないので、ggplotのgeom_histogramを使って ヒストグラムをプロットします。ほとんどRのggplot2と同じようにヒストプロットを表示する ことができます。
<ggplot: (8787705325877)> <ggplot: (8787705325877)> |
残念ながら、ggplotのヒストグラムにはバグがあるようで、最後の1列のプロットが正しく表示されません。
Saving 11.0 x 8.0 in image. Saving 11.0 x 8.0 in image. |
データのsummaryで平均値と分散が同じ値を示すことからデータ分布をポアソン分布と推定しています。
ポアソン分布では、データの平均がポアソン分布のλの値になることから、データの平均3.56のポアソン分布 を表示してみます。
yに0から9の値をセットし、r.dpoisでSageからRのdpoisを使って、Yに対するポアソン密度を計算し、その値(n())を浮動小数に変換して probにセットします。
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9] [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] |
[0.0284388247141845, 0.101242215982497, 0.180211144448844, 0.213850558079295, 0.190326996690573, 0.135512821643688, 0.0804042741752548, 0.0408913165805582, 0.0181966358783484, 0.00719778041410224] [0.0284388247141845, 0.101242215982497, 0.180211144448844, 0.213850558079295, 0.190326996690573, 0.135512821643688, 0.0804042741752548, 0.0408913165805582, 0.0181966358783484, 0.00719778041410224] |
Sageのlist_plotを使ってλ=3.56に対するポアソン分布をプロットします。
|
今後、ggplotを頻繁に活用して図化したいと考えているので、上記のポアソン分布をggplot のgeom_pointで表示してみます。
ggplotを使うことでとても簡単に高品質のグラフが出来上がります。
<ggplot: (8787705212321)> <ggplot: (8787705212321)> |
Saving 11.0 x 8.0 in image. Saving 11.0 x 8.0 in image. |
ggplotの大きな特徴の一つにデータの重ね合わせがあります。 ggplotのgeom_histogramとgeom_pointを足し合わせるだけで、グラフの重ね合わせができます。
<ggplot: (8787705160261)> <ggplot: (8787705160261)> |
Saving 11.0 x 8.0 in image. Saving 11.0 x 8.0 in image. |
ポアソン分布がλを変えることでどのように分布が変わるのか、見てみましょう。 SageのGraphisもグラフの重ね合わせたggplot同様にともて簡単です。
ですから、ggplotとSageの相性もとても良いのです。
|
|
尤度を「あてはまりの良さ」であり、L(λ)と書きます。ポアソン分布の尤度は、以下の式で与えられます。 $$ L(\lambda) = \prod_{i} p(y_i | \lambda) = \prod_{i} \frac{\lambda^{y_i} exp(-\lambda)}{y_i !} $$
L(λ) の対数とったものを対数尤度と呼び、上記の定義から以下のように表されます。 $$ log L(\lambda) = \sum_{i} \left( y_i log \lambda - \lambda - \sum_{k}^{y_i} log k \right) $$
このL(λ)が最大になるλは、L(λ)をλで偏微分することで求まります。 $$ \frac{\partial log L(\lambda)}{\partial \lambda} = \sum_i \left\{ \frac{y_i}{\lambda} - 1 \right \} = \frac{1}{\lambda} \sum_i y_i - k = 0 $$ この結果は、データの標本平均$\hat{\lambda}$となり、例題のデータに当てはめると以下の様になります。 $$ \hat{\lambda} = \frac{1}{50} \sum_i y_i = データの標本平均 = 3.56 $$
この式を例題のデータに対して、λ = 3.6としてその値の和をとると、図2.7のlogL = -97.3の値を得ます。
[1] -97.25516 [1] -97.25516 |
そこで、例題データに対する対数尤度を計算するlogLを以下のように定義し、 対数尤度が最大になるようなλを求めてみます。
|
確認のため、 今度は$\lambda=2.0の対数尤度は、logL = -121.9と一致します。
-121.881181787 -121.881181787 |
λの最尤推定値をグラフから求めてみます。λ=2.0から5.0を0.1刻みで計算してみます。
|
結果は、以下の様になり、標本の平均値3.56近辺で最大になっていることがわかります。
|
久保本のすごいところは、疑似乱数を使って最尤推定値のばらつきまでみているところです。
平均3.5のポアソン分布の50個のサンプルを1000セット作って、そのヒストグラムをプロット してみます。
|
今度は、Pandasの持つ、ヒストグラムプロット機能を使って最尤推定値のばらつきをプロットしてみます。
lambda 0 3.04 1 3.34 2 3.86 3 3.64 4 3.14 lambda 0 3.04 1 3.34 2 3.86 3 3.64 4 3.14 |
|
<ggplot: (8787703630761)> <ggplot: (8787703630761)> |
Saving 11.0 x 8.0 in image. Saving 11.0 x 8.0 in image. |
|