# # この企画は、雑誌や教科書にでているグラフをSageで再現し、 # グラフの意味を理解すると共にSageの使い方をマスターすることを目的としています。 #
## 前回に続き、データ解析のための統計モデリング入門 # （以下、久保本と書きます） # の第４章の例題をSageを使って再現してみます。 #
## 4章のメインは、久保本図4.9（以下に引用）のバイアス補正とその分布にあると思います。 #
#
#
#
# 数式処理システムSageのノートブックは、計算結果を表示するだけではなく、実際に動かすことができるの大きな特徴です。 # この機会にSageを分析に活用してみてはいかがでしょう。 #
# # # ## 最初に必要なライブラリーやパッケージをロードしておきます。jupyter用に新しくなったRUtil.pyも使います。 #
# # In[1]: get_ipython().run_cell_magic('HTML', '', '\n') # In[2]: # RとPandasのデータフレームを相互に変換する関数を読み込む # Rの必要なライブラリ r('library(ggplot2)') r('library(jsonlite)') # python用のパッケージ import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns get_ipython().run_line_magic('matplotlib', 'inline') # jupyter用のdisplayメソッド from IPython.display import display, Latex, HTML, Math, JSON # sageユーティリティ load('script/sage_util.py') # Rユーティリティ load('script/RUtil.py') # ## 「あてはまりの良さとモデルの複雑さ」を説明するために、３章のデータを使ってk=1とk=7の結果の比較が示されています。 # このような図を簡単に再現できれば、実際の分析の時にも役立つと思います。 #
## 3章のデータをも一度読み込み、statsmodelsのglmを使ってk=1とその結果をプロットします。 #
# # In[3]: # 3章のデータをもう一度読み込む d = pd.read_csv('http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/poisson/data3a.csv') d_C = d[d['f'] == 'C'] d_T = d[d['f'] == 'T'] # In[4]: # statsmodelsを使ってglmを計算します import statsmodels.formula.api as smf import statsmodels.api as sm from scipy.stats.stats import pearsonr # In[5]: # k=1でのglm回帰を実行 fit_1 = smf.glm('y ~ 1', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() # λの予測値をλ_1にセット d['lam_1'] = fit_1.predict() # In[6]: # λの予測値と一緒にプロット ax = d.sort_values(by=['x']).plot(x='x', y='lam_1') d_T.plot(kind='scatter', x='x', y='y', c='r', ax=ax, label='T') d_C.plot(kind='scatter', x='x', y='y', c='b', ax=ax, label='C') plt.legend(loc='upper left') plt.show() # ## 次にk=7（$log \lambda = \beta_1 + \beta_2 x + ... + \beta_7 x^6$）の計算は、$x^n$ # を計算するカスタム関数nPwOfをformulaの中で使用することで、実現できました。 #
## 予測値は、データのある点のみを結線したため、カクカクしているところもありますが、データの中心部分を通っています。 #
## k=1とk=7の結果をみるとk=7の結果は、あまりにも複雑な曲線を示しており、機械学習でいうところの「過学習」 # （Overfitting）の状態になっています。過学習では、学習用データ（ここでいうサンプルデータ）に引きずられ、 # 真のモデルから離れてしまいます。 #
# # In[7]: # カスタム関数nPwOfを使って多項式近似を行う def nPwOf(n, x): return x^n formula = 'y ~ nPwOf(7, x) + nPwOf(6, x) + nPwOf(5, x) + nPwOf(4, x) + nPwOf(3, x) + nPwOf(2, x) + x' fit_1 = smf.glm(formula=formula, data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() # λの予測値をλ_1にセット d['lam_7'] = fit_1.predict() # In[8]: # λの予測値と一緒にプロット ax = d.sort_values(by=['x']).plot(x='x', y='lam_7') d_T.plot(kind='scatter', x='x', y='y', c='r', ax=ax, label='T') d_C.plot(kind='scatter', x='x', y='y', c='b', ax=ax, label='C') plt.legend(loc='upper left') plt.show() # ## 逸脱度は、「あてはまりの悪さ」を表現する指標と久保本では説明されています。 #
## 最大対数尤度を$log L^*$とすると、逸脱度(D)は以下の様に定義されます。 # $$ # D = -2 log L^* # $$ #
## 4章で登場する逸脱度を以下の表に示します。 # 最大の逸脱度は、一定モデルで、最小の逸脱度はフルモデルです。 #
# # In[9]: # 様々な逸脱度 hdr = ['名前', '定義'] tbl = [ ['逸脱度(D)', '-2 log L*'], ['最小の逸脱度', 'フルモデルをあてはめたときのD'], ['残差逸脱度', 'D-最小のD'], ['最大の逸脱度', 'NullモデルをあてはめたときのD'], ['Null逸脱度', '最大のD - 最小のD']] Table2Html(tbl, hdr) # ## 各モデルのglmのfit結果を求めて、その結果を表にまとめます。 #
## ここで出てくるモデル選択の基準となるAIC(Akaike's information criterion)は、以下の様に定義されています。 # $$ # AIC = -2 { (最大対数尤度) - (最尤推定したパラメータの数)} = D + 2k # $$ # #
# # In[10]: # 逸脱度（Deviance）は、-2 log L* と定義されており、glmの結果の-2*llfから計算できる # 各モデルに対するfitを計算してDevianceの違いを整理してみる fit_f = smf.glm('y ~ f', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() fit_x = smf.glm('y ~ x', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() fit_x_f = smf.glm('y ~ x + f', data=d, family=sm.families.Poisson(link=sm.families.links.log)).fit() # ## データの個数と同じ数のパラメータを使ってあてはめたモデルがフルモデルです。 # これは、実際にそのようなモデルを作って計算するのでは無く、フルモデルはすべての観測値が予測値と一致するモデルであり、 # その結果、残差は０であり、対数尤度は、以下のようになります。 #
# # In[11]: # フルモデルの対数尤度の計算 fullLogL = sageobj(sum([r.dpois(y, y, log=True) for y in d.y])); fullLogL # ## 結果をみやすくするために、小数点以下3桁を表示する関数prf3を定義します。 #
# # In[12]: # 結果を小数点3桁で表示する def prf3(f): return '%.3f'%f # In[13]: # フルモデルは、k=Nで結果がすべて観測値と一致します。従って残差は0となります。 # log L*: fit.llf, residual deviance: fit.deviance, AIC: fit.aicで取得できる hdr = ['モデル', 'k', "log L*", 'deviance -2 log L*', 'residual deviance', 'AIC'] tbl = [['一定', '1', prf3(fit_1.llf), prf3(-2*fit_1.llf), prf3(fit_1.deviance), prf3(fit_1.aic)], ['f', '2', prf3(fit_f.llf), prf3(-2*fit_f.llf), prf3(fit_f.deviance), prf3(fit_f.aic)], ['x', '2', prf3(fit_x.llf), prf3(-2*fit_x.llf), prf3(fit_x.deviance), prf3(fit_x.aic)], ['x+f', '3', prf3(fit_x_f.llf), prf3(-2*fit_x_f.llf), prf3(fit_x_f.deviance), prf3(fit_x_f.aic)], ['フル', '100', prf3(fullLogL), prf3(-2*fullLogL), '0', prf3(-2*fullLogL+2*100)]] Table2Html(tbl, hdr) # ## いよいよ図４．９の計算に入ります。久保本ではモデル選択の指標としてAICが使えるのかを一定モデルを使って # 数値例で示しています。（これが久保本のすごいところ！） #
## 平均の対数尤度は、真のモデルを使って生成されたサンプルから求めた対数尤度の平均値であり、実際のモデルづくりでは # このようなことは不可能ですが、平均対数尤度と推定されたモデルの関係とこの差（バイアス）bの分布からAICの意味するところを # 表現しています。bの定義は以下の通りです。 # $$ # b = log L^* - E(Log L) # $$ #
# # ## 図４．９の(A)観測データが一つの場合について、真のモデルλ=8のポアソン分布から50個のサンプルを生成します。 #
## 準備としてサンプル生成関数genSample, 対数尤度を計算するlogLを定義します。 #
# # In[14]: # サンプル生成λ=8ポアソン分布のサンプルをN=50個生成する # sample = [ float(y.n()) for y in r.rpois(N, lambda_true)] # sageobjを使って変換する方が速いのでgenSample関数を使用する def genSample(lambda_target, sample_size): return sageobj(r('rpois(%d, %f)'%(sample_size, lambda_target))) # In[15]: # Rの乱数の種をセット r('set.seed(101)') # 真のλ=8から50個のポアソン分布のサンプルを生成する N = 50 lambda_true = 8 sample = genSample(lambda_true, N) sampleDf = pd.DataFrame({'y' : np.array(sample)}) # In[16]: # 対数尤度の計算 def logL(sample, lambda_estimated): return sageobj(sum([r.dpois(y, lambda_estimated, log=True) for y in sample])) # ## glmを使って一定モデルの結果を取得します。ここでは対数尤度が-116.7、 $\beta_1$の値は2.06と求まりました。 #
# # In[17]: # k=1でのglm回帰を実行 fit = smf.glm('y ~ 1', data=sampleDf, family=sm.families.Poisson(link=sm.families.links.log)).fit() # In[18]: # fit.summary2() # 最大対数尤度 mxLogL = fit.llf # 切片 beta_1 = fit.params[0] print mxLogL, beta_1 # ## この結果を図化し、logL_pltにセットします。 #
# # In[19]: # β= [1.95, 2.20]の範囲の対数尤度を求める logL_lst = [ (x,logL(sample, exp(x)) )for x in np.arange(1.95, 2.20, 0.01)] # log L*の曲線プロット logL_plt = list_plot(logL_lst, figsize=4, plotjoined =True) # ## つぎに、真のモデルから50個のサンプルを200セット生成し、平均の対数尤度を計算します。 #
# # In[20]: # 真のλ=8から生成されたサンプル200セットにβの推定値を使って対数尤度を計算 # logL200 = [ logL(genSample(lambda_true, N), exp(beta_1)) for i in range(50)] # とても遅いので、 logL200 = [ sageobj(r('sum(log(dpois(rpois(%d, %f), %f)))'%(N, lambda_true, exp(beta_1)))) for i in range(200)] # ## 求まった結果を一つにプロットします。logL*を赤○で、平均対数尤度を緑○で、各対数尤度を小さな青○で示し、 # 各βの対数尤度曲線と一緒に表示したのが、以下の図です。 #
# # In[21]: sample_plt = list_plot([ (beta_1, y) for y in logL200]) logL_pt = point((beta_1, mxLogL), size=40, color="red") logL_mean_pt = point((beta_1, mean(logL200)), size=40, color="green") (logL_plt + sample_plt + logL_pt + logL_mean_pt).show() # In[22]: # 実際の解析では不可能ですが、今回は真のモデルが分かっているので、評価用のデータを生成し、E(logL)を計算できます # ## 次に上記の処理を繰り返した図４．９(B)を試してみます。 # Sageでは、一度試した処理をコピー＆ペーストしてループに入れたり、関数にすることで結果を確認しながら # 処理を進めることで効率良く目的のグラフを得ることができます。 #
## 最初に空のGraphics=gを生成し、上記の処理と同じことをして、プロット結果をgに加えていきます。 # たったこれだけで、図４．９の(B)が再現できます。 #
## Rのseedをセットすることで、久保本の平均対数尤度に近い曲線が表示されました。（2016/07/17） #
# # In[23]: # 上記の計算を25回繰り返しプロット図4.9(B) g = Graphics() for i in range(25): sample = genSample(lambda_true, N) sampleDf = pd.DataFrame({'y' : np.array(sample)}) # k=1でのglm回帰を実行 fit = smf.glm('y ~ 1', data=sampleDf, family=sm.families.Poisson(link=sm.families.links.log)).fit() # 最大対数尤度 mxLogL = fit.llf # 切片 beta_1 = fit.params[0] logL200 = [ sageobj(r('sum(log(dpois(rpois(%d, %f), %f)))'%(N, lambda_true, exp(beta_1)))) for i in range(200)] # プロット g += list_plot([ (beta_1, y) for y in logL200]) g += point((beta_1, mxLogL), size=40, color="red") g += point((beta_1, mean(logL200)), size=40, color="green") # 結果を表示 g.show(figsize=4) # ## 非常に時間がかかったので試しに10セットのデータに対してバイアスbを計算し、その分布を表示してみました。 #
## バイアスbは、2.7と久保本の1.01とはずれています。 #
# # In[24]: # バイアスbの計算 def bias(lambda_true, N): sample = genSample(lambda_true, N) sampleDf = pd.DataFrame({'y' : np.array(sample)}) # k=1でのglm回帰を実行 fit = smf.glm('y ~ 1', data=sampleDf, family=sm.families.Poisson(link=sm.families.links.log)).fit() # 最大対数尤度 mxLogL = fit.llf # 切片 beta_1 = fit.params[0] logL200 = [ sageobj(r('sum(log(dpois(rpois(%d, %f), %f)))'%(N, lambda_true, exp(beta_1)))) for i in range(200)] return (mxLogL - sageobj(mean(logL200)) ) # In[25]: # 時間がかかるので、50個だけサンプルを生成しました bias_sample = [bias(lambda_true, N) for i in range(50)] # In[26]: mean(bias_sample) # In[27]: # biasのヒストグラムを表示する histogram(bias_sample, figsize=4) # ## 久保本のサポートページにあるRの関数を使ってサンプル1000セットでバイアスbを計算してみました。 # これでも久保本の1.01とはずれており、その分布をみると-20から20の区間でかなりぶれていることがわかります。 #
## バイアスの定義を変形すると、平均対数尤度$E(logL)$は以下の様になります。 # $$ # E(log L) = log L^* - b # $$ # ここで、一定モデルのk=1であることからバイアスbとkを入れ替えるとAICの定義となり、 # AICが予測の「悪さ」を表す指標となっていることがうなずけます。 # $$ # AIC = -2 \times (log L^* -1) # $$ #
# # # Rの関数をdata/bias.txtに定義します。 # In[28]: get_ipython().run_cell_magic('writefile', 'data/bias.txt', 'bias <- function(lambda.true, sample.size){\n sample.rpois <- rpois(sample.size, lambda.true)\n fit <- glm(sample.rpois~1, family=poisson) \n lambda.estimated <- exp(coef(fit))\n likelihood.mean <- mean(sapply(1:200, function(i){sum(log(dpois(rpois(N, lambda.true), lambda.estimated)))}))\n logLik(fit) - likelihood.mean\n}\n') # In[29]: # Rでbiasを計算する関数を定義すると非常に速く計算できる r('source("data/bias.txt") ') # サンプリングが1000くらいでもbiasの平均は、結構ぶれる r('N <- 50') r('lambda.true <- 8') r('bias.sampled <- sapply(1:1000, function(i) bias(lambda.true, N))') r('mean(bias.sampled)') # In[30]: graph = preGraph("images/fig-4.10-R.pdf") r('p <- qplot(bias.sampled, geom = "blank") + geom_histogram(aes(y = ..density..), colour = "black", fill = "white") + geom_density(alpha = 0.2, fill = "#6666FF") + geom_vline(aes(xintercept = mean(bias.sampled)), color = "red", linetype = "dashed", size = 2) ') r('plot(p)') postGraph(graph) # In[ ]: