[[FrontPage]]
#contents
2014/03/16からのアクセス回数 &counter;
ここで紹介したSageワークシートは、以下のURLからダウンロードできます。
http://www15191ue.sakura.ne.jp:8000/home/pub/38/
また、Sageのサーバを公開しているサイト（http://www.sagenb.org/, http://www15191ue.sakura.ne.jp:8000/）にユーザIDを作成することで、ダウンロードしたワークシートを
アップロードし、実行したり、変更していろいろ動きを試すことができます。
* Sageでグラフを再現してみよう：データ解析のための統計モデリング入門第２章 [#w4b6a337]
この企画は、雑誌や教科書にでているグラフをSageで再現し、
グラフの意味を理解すると共にSageの使い方をマスターすることを目的としています。
今回は、
[[データ解析のための統計モデリング入門>http://www.amazon.co.jp/dp/400006973X/]]
（以下、久保本と書きます）
の第２章の例題をSageを使って再現してみます。久保氏のブログは、WinBUGSの使い方などでよく拝見していましたが、
「データ解析のための統計モデリング入門」は自然科学を学ぶすべての人に薦めたい一冊です。まさに目から鱗の固まりです。
数式処理システムSageのノートブックは、計算結果を表示するだけではなく、実際に動かすことができるの大きな特徴です。
この機会にSageを分析に活用してみてはいかがでしょう。
** 前準備 [#c93b62ac]
最初に必要なライブラリーやパッケージをロードしておきます。新しくなったRUtil.pyも使います。
sageへの入力:
#pre{{
# 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')
}}
** 例題のデータ [#e9cad38b]
2章の例題に使われているデータは、架空の植物の第i番目の個体の種子数$y_i$を扱っています。
本の図２．１を以下に引用します。
&ref(Fig2.1.png);
2章のデータは、ネット公開されており、
[[本のサポートサイト>http://hosho.ees.hokudai.ac.jp/~kubo/ce/IwanamiBook.html]]
からダウンロードできます。
ここでは、DATAディレクトリにdata.RDataをアップし、このノートブックで参照できるようにしました。
データをロードしたら、summary関数とvar関数を使ってデータの素性を確認します。
ここで、確認することは、以下の3項目です。
- データの個数と欠損値の有無
- データの平均
- データの分散
sageへの入力:
#pre{{
# 2章のデータdata.RDataをダウンロードし、Rにロードする
r('load("%s")' % (DATA + 'data.RData'))
# summaryで平均と分散が大体同じであることがポアソン分布の特徴
print r('summary(data)')
print r('var(data)')
}}
#pre{{
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.00 2.00 3.00 3.56 4.75 7.00
[1] 2.986122
}}
*** データの可視化 [#y496e6c7]
データの分布を見るには、Rのtable関数を使って度数分布を計算すると便利です。またグラフに表示するときには、hist関数を使用します。
sageへの入力:
#pre{{
# 度数分布を計算
r('table(data)')
}}
#pre{{
data
0 1 2 3 4 5 6 7
1 3 11 12 10 5 4 4
}}
sageへの入力:
#pre{{
# Rのhistを使ってヒストグラムをプロットする
graph = preGraph("fig2.2-R.png")
r('p <- hist(data, breaks = seq(-0.5, 9.5, 1))')
r('plot(p)')
postGraph(graph)
}}
&ref(fig2.2-R.png);
** Pandasを使って同様の処理に挑戦 [#ab6eb20c]
Rのライブラリはとても豊富で実績もありますが、個人的には関数の定義が分かりずらいのであまり好きではありません。
そこで、SageにPandasライブラリをインストールして、これを使って同じような処理をしてみました。
Rから配列データを取り出すには、sageobj関数を使います。
sageへの入力:
#pre{{
# 同様の処理をPandasを使ってやってみます
# rからSageにデータを変換
data = sageobj(r('data')); data
}}
#pre{{
[2, 2, 4, 6, 4, 5, 2, 3, 1, 2, 0, 4, 3, 3, 3, 3, 4, 2, 7, 2, 4, 3, 3, 3, 4, 3, 7, 5, 3, 1, 7, 6, 4, 6, 5, 2, 4, 7, 2, 2, 6, 2, 4, 5, 4, 5, 1, 3, 2, 3]
}}
*** データフレームに変換 [#xbd16369]
Rから取り込んだ2章のデータをPandasのデータフレームにします。
sageへの入力:
#pre{{
# ggplotでプロットできるようにDataFrameにする
orgData = pd.DataFrame(data, columns = ['data'])
}}
データフレームにしたorgDataからRのsummaryと同様の情報をdescribeを使って得ることができます。
sageへの入力:
#pre{{
# RのsummaryのようにPandasのDataFrameの情報を出力するには、describeを使います
orgData.describe()
}}
#pre{{
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
[8 rows x 1 columns]
}}
*** ggplotを使ってヒストグラムを表示 [#zd59b702]
Sageにはヒストグラムをプロットする関数がないので、ggplotのgeom_histogramを使って
ヒストグラムをプロットします。ほとんどRのggplot2と同じようにヒストプロットを表示する
ことができます。
sageへの入力:
#pre{{
# Sageにはヒストグラムをプロットする関数がないので、ggplotを使用する
# 50個の度数になるようにprob*50としてプロットしている
ggplot(orgData, aes(x='data')) + geom_histogram(binwidth=1)
}}
#pre{{
<ggplot: (15907857)>
}}
残念ながら、ggplotのヒストグラムにはバグがあるようで、最後の1列のプロットが正しく表示されません。
sageへの入力:
#pre{{
# ヒストグラムの最後のプロットがRと異なる
ggsave('fig_2.2-g.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig_2.2-g.png);
** ポアソン分布の表示 [#z55d7f88]
データのsummaryで平均値と分散が同じ値を示すことからデータ分布をポアソン分布と推定しています。
ポアソン分布では、データの平均がポアソン分布のλの値になることから、データの平均3.56のポアソン分布
を表示してみます。
yに０から9の値をセットし、r.dpoisでSageからRのdpoisを使って、Yに対するポアソン密度を計算し、その値（n()）を浮動小数に変換して
probにセットします。
sageへの入力:
#pre{{
# ポアソン分布を表示
y = range(10); y
}}
#pre{{
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
}}
sageへの入力:
#pre{{
# dpoisの戻り値が複素数になっているので、floatで浮動小数に変換
prob = [float(r.dpois(v, 3.56).n()) for v in y]; prob
}}
#pre{{
[0.0284388247141845, 0.101242215982497, 0.180211144448844, 0.213850558079295, 0.190326996690573,
0.135512821643688, 0.0804042741752548, 0.0408913165805582, 0.0181966358783484, 0.00719778041410224]
}}
*** ポアソン分布をプロットする [#s2f9884a]
Sageのlist_plotを使ってλ=3.56に対するポアソン分布をプロットします。
sageへの入力:
#pre{{
# Sage版 Fig. 2.4
list_plot(zip(y, prob), figsize=5)
}}
&ref(sage1.png);
*** 同様の図をggplotで表示してみる [#s77ae509]
今後、ggplotを頻繁に活用して図化したいと考えているので、上記のポアソン分布をggplot
のgeom_pointで表示してみます。
ggplotを使うことでとても簡単に高品質のグラフが出来上がります。
sageへの入力:
#pre{{
# yとprobからデータフレームを作る
# 同様のプロットをggplotで表示してみる
df = pd.DataFrame(zip(y, prob), columns = ['y', 'prob'])
# ggplotでプロットしてみる
ggplot(df, aes(x='y', y='prob')) + geom_point()
}}
#pre{{
<ggplot: (15979209)>
}}
sageへの入力:
#pre{{
ggsave('fig_2.4.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig_2.4.png);
*** 観測データとの重ね合わせ [#rc2c3349]
ggplotの大きな特徴の一つにデータの重ね合わせがあります。
ggplotのgeom_histogramとgeom_pointを足し合わせるだけで、グラフの重ね合わせができます。
sageへの入力:
#pre{{
# 観測データと平均3.56のポアソン分布を重ね合わせてみる
# 50個の度数になるようにprob*50としてプロットしている
ggplot(orgData, aes(x='data')) + \
geom_histogram(binwidth=1, alpha=0.6, color='grey') + \
geom_point(df, aes(x='y', y='prob*50'), size=100, color="darkred")
}}
#pre{{
<ggplot: (16444061)>
}}
sageへの入力:
#pre{{
ggsave('fig_2.5.png', dpi=50)
}}
#pre{{
Saving 11.0 x 8.0 in image.
}}
&ref(fig_2.5.png);
** ポアソン分布の特徴 [#n22f16f6]
ポアソン分布がλを変えることでどのように分布が変わるのか、見てみましょう。
SageのGraphisもグラフの重ね合わせたggplot同様にともて簡単です。
ですから、ggplotとSageの相性もとても良いのです。
sageへの入力:
#pre{{
# ポアソン関数p(y, λ)を定義
def p(y, lam):
return lam^y*e^(-lam)/factorial(y)
}}
sageへの入力:
#pre{{
# 様々な平均（λ）のポアソン分布
# 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)
}}
&ref(sage2.png);
** ポアソン分布の最尤推定 [#oabe9ecc]
尤度を「あてはまりの良さ」であり、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の値を得ます。
sageへの入力:
#pre{{
# lambda=3.6の対数尤度を計算してみる。本の図2.7のlog Lと同じ値であることを確認
orgData.data.apply(lambda y : r.dpois(y, 3.6, log=True)).sum()
}}
#pre{{
[1] -97.25516
}}
*** 例題データの対数尤度を計算する関数logLを定義 [#h6ff8e9b]
そこで、例題データに対する対数尤度を計算するlogLを以下のように定義し、
対数尤度が最大になるようなλを求めてみます。
sageへの入力:
#pre{{
# 対数尤度の計算
def logL(m):
return orgData.data.apply(lambda y : r.dpois(y, m, log=True)).sum().n()
}}
確認のため、 今度は$\lambda=2.0の対数尤度は、logL = -121.9と一致します。
sageへの入力:
#pre{{
# 関数の確認、今度はlambda=2の値でチェック
print logL(2)
}}
#pre{{
-121.881181786917
}}
*** 例題の最尤推定値 [#he30177e]
λの最尤推定値をグラフから求めてみます。λ=2.0から5.0を0.1刻みで計算してみます。
sageへの入力:
#pre{{
# グラフから最大となる対数尤度を求める（少し時間が掛かる）
lams = np.arange(2.0, 5.0, 0.1)
Ls = [logL(x) for x in lams]
}}
結果は、以下の様になり、標本の平均値3.56近辺で最大になっていることがわかります。
sageへの入力:
#pre{{
# plotだと時間がかかるので、list_plotで代用
# Fig. 2.8
list_plot(zip(lams, Ls), figsize=5, plotjoined =True)
}}
&ref(sage3.png);
** 疑似乱数と最尤推定値のばらつき [#uae72a6b]
久保本のすごいところは、疑似乱数を使って最尤推定値のばらつきまでみているところです。
平均3.5のポアソン分布の50個のサンプルを1000セット作って、そのヒストグラムをプロット
してみます。
sageへの入力:
#pre{{
# 平均3.5のポアソン分布のサンプルを50個を1000セット生成して、最尤推定値のばらつきを見る
poisSet3000 = [ r.rpois(50, 3.5).mean().n() for i in range(1000)]
}}
*** 最尤推定値のばらつき [#fab7402c]
今度は、Pandasの持つ、ヒストグラムプロット機能を使って最尤推定値のばらつきをプロットしてみます。
sageへの入力:
#pre{{
# Pandasのグラフ機能を使ってヒストグラムを表示する（ggplotのヒストグラムはバグありみたいだから）
df3000 = pd.DataFrame({ 'lambda': np.array(poisSet3000)}); df3000.head()
}}
#pre{{
lambda
0 (3.7+0j)
1 (3.36+0j)
2 (3.58+0j)
3 (3.78+0j)
4 (3.56+0j)
[5 rows x 1 columns]
}}
sageへの入力:
#pre{{
df3000.hist()
plt.savefig("df3000.png", dpi=70)
}}
&ref(df3000.png);
** コメント [#p46ea21d]
#vote(おもしろかった[2],そうでもない[0],わかりずらい[0])
#vote(おもしろかった[2],そうでもない[1],わかりずらい[0])
皆様のご意見、ご希望をお待ちしております。
- 「データ解析のための統計モデリング入門」の著者の久保拓哉氏にリツイートして頂いてとても嬉しいです。 -- [[竹本 浩]] &new{2014-03-17 (月) 12:27:48};
- 偏微分が０であることが式に追加した -- [[竹本 浩]] &new{2014-03-26 (水) 16:37:15};
#comment_kcaptcha
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
