Hiroshi TAKEMOTO (take@pwv.co.jp)

Sageでグラフを再現してみよう：データ解析のための統計モデリング入門第２章

この企画は、雑誌や教科書にでているグラフをSageで再現し、 グラフの意味を理解すると共にSageの使い方をマスターすることを目的としています。

今回は、データ解析のための統計モデリング入門 （以下、久保本と書きます） の第２章の例題をSageを使って再現してみます。久保氏のブログは、WinBUGSの使い方などでよく拝見していましたが、 「データ解析のための統計モデリング入門」は自然科学を学ぶすべての人に薦めたい一冊です。まさに目から鱗の固まりです。

数式処理システムSageのノートブックは、計算結果を表示するだけではなく、実際に動かすことができるの大きな特徴です。 この機会にSageを分析に活用してみてはいかがでしょう。

前準備

最初に必要なライブラリーやパッケージをロードしておきます。新しくなったRUtil.pyも使います。

{{{id=1| # Rの必要なライブラリ #r("install.packages('ggplot2')") r('library(ggplot2)') r('library(jsonlite)') # python用のパッケージ import pandas as pd import numpy as np import matplotlib.pyplot as plt from ggplot import * # RUtilにRとPandasのデータフレームを相互に変換する関数を追加 load(DATA + 'RUtil.py') /// }}}

例題のデータ

2章の例題に使われているデータは、架空の植物の第i番目の個体の種子数$y_i$を扱っています。 本の図２．１を以下に引用します。

2章のデータは、ネット公開されており、本のサポートサイトからダウンロードできます。 ここでは、DATAディレクトリにdata.RDataをアップし、このノートブックで参照できるようにしました。

データをロードしたら、summary関数とvar関数を使ってデータの素性を確認します。 ここで、確認することは、以下の3項目です。

{{{id=62| # 2章のデータdata.RDataをダウンロードし、Rにロードする r('load("%s")' % (DATA + 'data.RData')) # summaryで平均と分散が大体同じであることがポアソン分布の特徴 print r('summary(data)') print r('var(data)') /// 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関数を使用します。

{{{id=63| # 度数分布を計算 r('table(data)') /// data 0 1 2 3 4 5 6 7 1 3 11 12 10 5 4 4 }}} {{{id=64| # Rのhistを使ってヒストグラムをプロットする graph = preGraph("fig2.2-R.pdf") r('p <- hist(data, breaks = seq(-0.5, 9.5, 1))') #r('plot(p)') postGraph(graph) /// }}}

Pandasを使って同様の処理に挑戦

Rのライブラリはとても豊富で実績もありますが、個人的には関数の定義が分かりずらいのであまり好きではありません。

そこで、SageにPandasライブラリをインストールして、これを使って同じような処理をしてみました。

Rから配列データを取り出すには、sageobj関数を使います。

{{{id=17| # 同様の処理をPandasを使ってやってみます # rからSageにデータを変換 data = sageobj(r('data')) /// }}}

データフレームに変換

Rから取り込んだ2章のデータをPandasのデータフレームにします。

{{{id=18| # ggplotでプロットできるようにDataFrameにする orgData = pd.DataFrame(data, columns = ['data']) /// }}}

データフレームにしたorgDataからRのsummaryと同様の情報をdescribeを使って得ることができます。

{{{id=60| # RのsummaryのようにPandasのDataFrameの情報を出力するには、describeを使います orgData.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 }}}

ggplotを使ってヒストグラムを表示

Sageにはヒストグラムをプロットする関数がないので、ggplotのgeom_histogramを使って ヒストグラムをプロットします。ほとんどRのggplot2と同じようにヒストプロットを表示する ことができます。

{{{id=19| # Sageにはヒストグラムをプロットする関数がないので、ggplotを使用する # 50個の度数になるようにprob*50としてプロットしている ggplot(orgData, aes(x='data')) + geom_histogram(binwidth=1.0) /// }}}

残念ながら、ggplotのヒストグラムにはバグがあるようで、最後の1列のプロットが正しく表示されません。

{{{id=20| # ヒストグラムの最後のプロットがRと異なる GgSave('fig_2.2-g.png') /// Saving 11.0 x 8.0 in image. }}}

ポアソン分布の表示

データのsummaryで平均値と分散が同じ値を示すことからデータ分布をポアソン分布と推定しています。

ポアソン分布では、データの平均がポアソン分布のλの値になることから、データの平均3.56のポアソン分布 を表示してみます。

yに０から9の値をセットし、r.dpoisでSageからRのdpoisを使って、Yに対するポアソン密度を計算し、その値（n()）を浮動小数に変換して probにセットします。

{{{id=4| # ポアソン分布を表示 y = range(10); y /// [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] }}} {{{id=5| # dpoisの戻り値が複素数になっているので、floatで浮動小数に変換 prob = [float(sageobj(r.dpois(v, 3.56))) for v in y]; prob /// [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に対するポアソン分布をプロットします。

{{{id=7| # Sage版 Fig. 2.4 list_plot(zip(y, prob), figsize=5) /// }}}

同様の図をggplotで表示してみる

今後、ggplotを頻繁に活用して図化したいと考えているので、上記のポアソン分布をggplot のgeom_pointで表示してみます。

ggplotを使うことでとても簡単に高品質のグラフが出来上がります。

{{{id=10| # yとprobからデータフレームを作る # 同様のプロットをggplotで表示してみる df = pd.DataFrame(zip(y, prob), columns = ['y', 'prob']) # ggplotでプロットしてみる ggplot(df, aes(x='y', y='prob')) + geom_point() /// }}} {{{id=11| GgSave('fig_2.4.png') /// Saving 11.0 x 8.0 in image. }}}

観測データとの重ね合わせ

ggplotの大きな特徴の一つにデータの重ね合わせがあります。 ggplotのgeom_histogramとgeom_pointを足し合わせるだけで、グラフの重ね合わせができます。

{{{id=65| # 観測データと平均3.56のポアソン分布を重ね合わせてみる # 50個の度数になるようにprob*50としてプロットしている df.prob = df.prob*50 ggplot(orgData, aes(x='data')) + \ geom_histogram(binwidth=1.0, alpha=0.6, color='grey') + \ geom_point(df, aes(x='y', y='prob'), size=100, color="darkred") /// }}} {{{id=66| GgSave('fig_2.5.png') /// Saving 11.0 x 8.0 in image. }}}

ポアソン分布の特徴

ポアソン分布がλを変えることでどのように分布が変わるのか、見てみましょう。 SageのGraphisもグラフの重ね合わせたggplot同様にともて簡単です。

ですから、ggplotとSageの相性もとても良いのです。

{{{id=37| # ポアソン関数p(y, λ)を定義 def p(y, lam): return lam^y*e^(-lam)/factorial(y) /// }}} {{{id=40| # 様々な平均（λ）のポアソン分布 # Fig. 2.6 g = Graphics() for lam, cls in [(3.5, 'red'), (7.7, 'blue'), (15.1, 'green')]: g = g + list_plot([p(y, lam) for y in range(20)], color=cls) g.show(figsize=5) /// }}}

ポアソン分布の最尤推定

尤度を「あてはまりの良さ」であり、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の値を得ます。

{{{id=43| # lambda=3.6の対数尤度を計算してみる。本の図2.7のlog Lと同じ値であることを確認 orgData.data.apply(lambda y : r.dpois(y, 3.6, log=True)).sum() /// [1] -97.25516 }}}

例題データの対数尤度を計算する関数logLを定義

そこで、例題データに対する対数尤度を計算するlogLを以下のように定義し、 対数尤度が最大になるようなλを求めてみます。

{{{id=44| # 対数尤度の計算 def logL(m): return sageobj(orgData.data.apply(lambda y : r.dpois(y, m, log=True)).sum()) /// }}}

確認のため、 今度は$\lambda=2.0の対数尤度は、logL = -121.9と一致します。

{{{id=47| # 関数の確認、今度はlambda=2の値でチェック print logL(2) /// -121.881181787 }}}

例題の最尤推定値

λの最尤推定値をグラフから求めてみます。λ=2.0から5.0を0.1刻みで計算してみます。

{{{id=50| # グラフから最大となる対数尤度を求める（少し時間が掛かる） lams = np.arange(2.0, 5.0, 0.1) Ls = [logL(x) for x in lams] /// }}}

結果は、以下の様になり、標本の平均値3.56近辺で最大になっていることがわかります。

{{{id=48| # plotだと時間がかかるので、list_plotで代用 # Fig. 2.8 list_plot(zip(lams, Ls), figsize=5, plotjoined =True) /// }}}

疑似乱数と最尤推定値のばらつき

久保本のすごいところは、疑似乱数を使って最尤推定値のばらつきまでみているところです。

平均3.5のポアソン分布の50個のサンプルを1000セット作って、そのヒストグラムをプロット してみます。

{{{id=51| # 平均3.5のポアソン分布のサンプルを50個を1000セット生成して、最尤推定値のばらつきを見る poisSet3000 = [ sageobj(r.rpois(50, 3.5).mean()) for i in range(1000)] /// }}}

最尤推定値のばらつき

今度は、Pandasの持つ、ヒストグラムプロット機能を使って最尤推定値のばらつきをプロットしてみます。

{{{id=57| # Pandasのグラフ機能を使ってヒストグラムを表示する df3000 = pd.DataFrame({ 'lambda': np.array(poisSet3000)}); df3000.head() /// lambda 0 3.04 1 3.34 2 3.86 3 3.64 4 3.14 }}} {{{id=58| df3000.hist() plt.savefig("df3000.png", dpi=70) /// }}} {{{id=67| # ggplotでもうまくでるようになりました（2015/08/09） # geom_histogramの例 ggplot(df3000, aes(x='lambda')) + geom_histogram(binwidth=0.2) /// }}} {{{id=70| GgSave("df3000-gg.png") /// Saving 11.0 x 8.0 in image. }}} {{{id=71| /// }}}