#!/usr/bin/env python
# coding: utf-8

# # NagoyaStat #2
# # データ解析のための統計モデリング入門　第４章

# ### GLMのモデル選択
# ### -AICとモデルの予測の良さ-
# + @nishiokya
# + 2016.10.01

# # 0.はじめに
# ##
# - モデルを選択する方法
# - AICを使い方を理解する
# - AICが正しい理由
# 
# ## 前提
# - AICの展開はしない

# # はじめに
# ## 自己紹介
# + @nishiokya
#     + エンジニア
#         + 統計検定２級(2014)
#         + AWS/Azure 
#         + 地図＋ビッグデータ
#     + Focus
#         + カーナビのセンサー情報の収集・変換・分析
# <div style="text-align: left;">
# <img src='https://2.bp.blogspot.com/_UBhXUUQQArM/S5IM58yOy8I/AAAAAAAAAbk/bN24ilKext0/S220-s73/787574_3864849165_bigger.jpg'>
# </div>

# # 0.はじめに
# ## 第４章の目次
# ### 4. GLMのモデル選択 −AICとモデルの予測の良さ-
# #### 4.1 データはひとつ，モデルはたくさん
# #### 4.2 統計モデルの当てはまりの悪さ：逸脱度
# #### 4.3 モデル選択基準AIC 
# #### 4.4 AICを説明するためのまた別の例題 
# #### 4.5 なぜAICでモデル選択をしてよいのか
# #### 4.6 この章のまとめと参考文献

# # 0.はじめに
# ## 振り返り
# ### ４章の立ち位置
# ![第2章以降の説明の流れ ](http://id.fnshr.info/wp-content/uploads/sites/2/2012/08/kubo_nagare.png "第2章以降の説明の流れ ")
# <div style="text-align: right;">
# まえがき P.viiより
# </div>

# # 0.はじめに
# ## 用語
# ### 第４章の用語集
# - AIC(赤池情報量規準 )
# - 逸脱度
#     - Null逸脱度
#     - 最大の逸脱度
#     - 最小の逸脱度
#     - 残差逸脱度
# - 平均対数尤度

# # ４章のサマリー
# - 4-2まで 当てはまりの良さと予測のよさ
#     - 3章で登場したGLM3種の統計モデルを比較する
#     - 逸脱度から統計モデルのパラメータ数と*あてはまり*の関係を理解する
# - 4-3 AICの使い方
#     - AICからパラメータ数と*予測* の関係を理解する
# - 4-4以降　AICの詳細
#     - AICがなぜモデル判定に利用できるかを理解する

# # ４章のサマリー
# ## 良いモデルとは
# ### *観測データ*にあてはまりが良いものが「良い」統計モデルか？
#  

# <p style="text-align:center;font-size:150%;line-height:150%;color:red">NO</p>
# 
# ###  母集団を「良く推測をするものが良いモデル」または、未知のデータを「良い予測をするものが良いモデル」
# ### 「良い予測をするものが良いモデル」という考えで、統計モデルを選択する方法をまなぶ
# 

# # ４章のサマリー
# ## 統計モデル選択
# ### 複数の統計モデルの中からなんらかの意味で「良い」モデルを選ぶこと
# - 良いモデルの選択基準は複数ある
#     - 良いあてはまりをする
#     - 良い予測をする

# # ４章のサマリー
# ## 良い予測とは
# + 同じデータ取得方法で異なるデータを得たとき、異なるデータを正確に言い当てれる

# # ４章のサマリー
# ## 良い予測するモデルを選ぶには
# ### 良い予測するをする選択基準を利用する
# - 良いモデルの判定方法は複数ある
#     - 「あてはまりのよい」モデルの選択基準　
#         - 最大尤度（最大対数尤度）
#         - p値（？）
#     - 「予測のよい」モデルの選択基準
#         - AIC

# # 4-1 データはひとつ，モデルはたくさん
# ## 3章で作成した統計モデルのうち一番「良い」モデルはどれ？
# <div style="text-align:center;font-size:150%;">3章の統計モデルの概要</div>
# 
# 
# モデル名  |  概要 | Rの構文|
# ---- | ---- | ----
#  xモデル  |  体サイズが影響するモデル | glm(y~x,data = d,family = poisson)
#  fモデル  |  施肥効果が影響するモデル | glm(y~f,data = d,family = poisson)
#  x+fモデル  |  体サイズと施肥効果が影響するモデル | glm(y~x+f,data=d,family=poisson)

# <div style="text-align:center;font-size:150%;">4章で追加される統計モデル</div>
# 
# モデル名  |  概要 | Rの構文|
# ---- | ----| ----
#   一定モデル*1  |  体サイズの効果も施肥効果も影響しないモデル。RではNULLモデルという| glm(y~1,data=d,family=poisson)
#  フルモデル  |  後述 |

# # 4-1 データはひとつ，モデルはたくさん
# ## 統計モデルを最大対数尤度で評価
# ### xモデルの最大対数尤度の算出方法
# 
# logLik(glm(y~x,data = d,family = poisson))
# 
# 
# <div style="text-align:center;font-size:120%;">モデルごとの最大対数尤度</div>
# 
# モデル名  |  最大対数尤度 |パラメータ数(k)
# ---- | ---- | ----
#  xモデル  |  -235.3863| 2
#  fモデル  |  -237.6273 | 2
#  x+fモデル  | -235.2937 | 3
#  一定モデル  | -237.6432  | 1
#  フルモデル  | -192.2  | 100

# ### モデルごとの最大対数尤度を比較する
# + x+fモデルの最大対数尤度が高い
#     + あてはまりは良い
# + パラメータ数が多い
#     + モデルが複雑
# <p style="text-align:center;font-size:150%;line-height:150%;color:red">この結果をどう評価するのか？</p>

# # 4.2 統計モデルのあてはまりの悪さ：逸脱度
# ## 逸脱度(deviance)とは
# ### 最大対数尤度を変形した統計量
# 
# ``` 
# D =  -2 logL* ・・・・　最大対数尤度に -2 をかけている
# ```
# 
# 

# # 4.2 統計モデルのあてはまりの悪さ：逸脱度
# ## 逸脱度とは
# ### 「あてはまりの良さ」ではなく「あてはまりの悪さ」で表現する指標
# 
# 
# <div style="text-align:center;font-size:150%;">モデルごとの逸脱度</div>
# 
# モデル名  |  最大対数尤度 |パラメータ数(k)| 逸脱度 | コメント
# ---- | ---- | ---- | ----| ----
#  xモデル  |  -235.3863| 2 |  470.7725 | 
#  fモデル  |  -237.6273 | 2 | 475.2546 |
#  x+fモデル  | -235.2937 | 3 | 470.5874 |
#  一定モデル  | -237.6432  | 1 | 475.2864 | 最大の逸脱度
#  フルモデル  | -192.2  | 100 | 384.4 |  最小の逸脱度
#  
# ### パラメータ数を増やせば逸脱度はさがる
# #### 　最大対数尤度に -2をかけただけなので当然

# <div style="text-align:center;font-size:180%;">対数尤度と逸脱度の関係</div>
# 
# 対数尤度(log likelihood)  |  逸脱度(deviance)
# ---- | ----
# あてはまりの良さ |  あてはまりの悪さ

# # 4.2 統計モデルのあてはまりの悪さ：逸脱度
# ## さまざまな逸脱度
# ### さまざまな逸脱度の定義
# 名前  |  定義 | glmのレスポンス
# ---- | ---- | ----
#  逸脱度 |当てはまりの悪さ(-2logL*)|
#  最小の逸脱度  |  フルモデルをあてはめたときの逸脱度 |
#  最大の逸脱度  |  一定モデルを当てはめたときの逸脱度 |
#  Null逸脱度  |  最大の逸脱度 - 最小逸脱度  | Null Deviance
#  残差逸脱度  |  逸脱度 - 最小逸脱度 |  Residual Deviance
#  logL*は最大対数尤度
#  

# # 4.2 統計モデルのあてはまりの悪さ：逸脱度
# ## 最小の逸脱度とは
# ### 「もっともあてはまりの良いモデル」の逸脱度
# ### 「フルモデル」が最小の逸脱度となる
# #### フルモデル
# #### データ１つ１つに対してパラメータを１つ１つ割り当てたモデル
# #### 推定用データに最適化しているため、予測には使えない

# # 4.2 統計モデルのあてはまりの悪さ：逸脱度
# ## 最大の逸脱度とは
# ### 「一定モデル（NULLモデル）」の逸脱度
# ### パラメータが１つしかない
# 

# # 4.2 統計モデルのあてはまりの悪さ：逸脱度
# ## Null逸脱度とは
# 最大の逸脱度（「一定モデル（NULLモデル）」の逸脱度） - 最小の逸脱度となる

# # 4.2 統計モデルのあてはまりの悪さ：逸脱度
# ## 残差逸脱度とは
# <div style="text-align:center;font-size:180%;">xモデルの逸脱度一覧</div>
# 
# 名前  |  定義 | 式
# ---- | ---- | ----
#  最大対数尤度 | -235.3863 |  logLik(glm(y~x,data = d,family = poisson))
#  逸脱度 | 470.7725  | -2 * logLik(glm(y~x,data = d,family = poisson))
#  残差逸脱度  | 84.99  |  470.7725 - 最小逸脱度
# 
# ```
# > glm(y~x,data=d,family=poisson)
# Call:  glm(formula = y ~ x, family = poisson, data = d)
# Coefficients:
# (Intercept)            x  
#     1.29172      0.07566  
# Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
# Null Deviance:	    89.51 
# Residual Deviance: 84.99 	AIC: 474.8
# ```

# # 4.2 統計モデルの当てはまりの悪さ：逸脱度
# ## 残差逸脱度とは
# <div style="text-align:center;font-size:180%;">モデルごとの残差逸脱度</div>
# 
# モデル名  |  最大対数尤度 |k| <p style="text-align:center;font-size:150%;line-height:150%;">逸脱度</p> |  <p style="text-align:center;font-size:150%;line-height:150%;">残差逸脱度</p>
# ---- | ---- | ---- | ----| ----| ----
#  NULLモデル(一定モデル）  | -237.6432  | 1 | <p style="text-align:center;font-size:150%;line-height:150%;">475.2864</p> | <p style="text-align:center;font-size:150%;line-height:150%;">89.5</p>
#  xモデル  |  -235.3863| 2 | <p style="text-align:center;font-size:150%;line-height:150%;"> 470.7725</p> |  <p style="text-align:center;font-size:150%;line-height:150%;">89.5 </p>
#  fモデル  |  -237.6273 | 2 |<p style="text-align:center;font-size:150%;line-height:150%;"> 475.2546</p> | <p style="text-align:center;font-size:150%;line-height:150%;">85.0</p>
#  x+fモデル  | <p style="color:red;">-235.2937</p> | 3 | <p style="text-align:center;font-size:150%;line-height:150%;color:red">470.5874 </p>| <p style="text-align:center;font-size:150%;line-height:150%;color:red">84.8</p>
#   フルモデル  | -192.2  | 100 | <p style="text-align:center;font-size:150%;line-height:150%;">384.4</p> | <p style="text-align:center;font-size:150%;line-height:150%;"> 0</p>
# 
# 

# # 4.2 統計モデルのあてはまりの悪さ：逸脱度
# ## 残差逸脱度は
# ### パラメータ数さえ増やせば残差逸脱度はどんどん小さくなりあてはまりが良くなる
# #### 残差逸脱度を０近くなるということはフルモデルにちかくなる？？
# ###  残差逸脱度ではパラメータ数が多いx+fモデル が良いモデルと選択

# # 4.3 モデル選択基準 AIC
# ## パラメータが多い方が、*当てはまりは良くなる*
# ### パラメータが多い方が、たまたま得られたデータが当てはまりが良くなる
# ## パラメータが多い方が、かならず*予測は良くなる*わけではない

# # 4.3 モデル選択基準 AIC
# ## モデル選択
# ### 複数の統計モデルの中からなんらかの意味で「良い」モデルを選ぶこと
# ### 「良い」モデルを選ぶ基準をモデル選択基準という
# ### AICは予測の良さを重視するモデル選択基準
# 

# # 4.3 モデル選択基準 AIC
# ## AICの計算式 
# ```
# AIC = -2 {(最大対数尤度) - 最尤推定したパラメータ数) }
#     = -2 ( logL* -K )
#     = D +2k
# ```
# 
# AICの値が一番小さいモデルが良い予測をするモデル
# 
# 

# # 4.3 モデル選択基準 AIC
# ## 統計モデルのAIC計算結果
# モデル名  |  最大対数尤度 |k| 逸脱度 | 残差逸脱度 | <p style="text-align:center;font-size:150%;line-height:120%;">AIC</p>
# ---- | ---- | ---- | ----| ----| ----
#  NULLモデル(一定モデル）  | -237.6432  | 1 | 475.2864 | 89.5 | <p style="text-align:center;font-size:120%;line-height:120%;">477.3<p>
#  xモデル  |  -235.3863| 2 |  470.7725 |  89.5 |<p style="text-align:center;font-size:150%;line-height:120%;color:red"> 474.8</p>
#  fモデル  |  -237.6273 | 2 | 475.2546 | 85.0| <p style="text-align:center;font-size:150%;line-height:120%;">479.3</p>
#  x+fモデル  | -235.2937 | 3 | 470.5874 | 84.8| <p style="text-align:center;font-size:150%;line-height:120%;">476.6</p>
#   フルモデル  | -192.2  | 50 | 384.4 |  0 | <p style="text-align:center;font-size:150%;line-height:120%;">585.8</p>
# ### AICでは xモデル をx+fモデルより良いモデルと選択

# # 4.4 AICを説明するための別の説明
# ## ネストしている関係のモデルとは
# ### 結論: ネストしている関係のモデルだとAICの精度はより良くなる
# 
# #### ネストした関係のモデル
# ```
# 1) logλi = β1 ・・・一定モデル(K=1)
# 2) logλi = β1 +  β2 * xi ・・・x モデル(K=2)
# ```
# + x モデルのβ2を0とすると一定モデルとおなじ。
#     + 1)と2) はネストしている関係のモデルという　

# # 4.5 なぜAICでモデル選択をして良いのか
# ## 4.5.1 統計モデルの予測の良さ
# ## 4.5.2 最大対数尤度のバイアス補正
# ## 4.5.3 ネストしているGLM間のAIC比較

# # 4.5.1 統計モデルの予測の良さ
# ## 平均対数尤度とは
# + 統計モデルに対して構築時と同じデータ取得方法で異なるデータセットを複数回あてはめた時の最大対数尤度の平均
#     + 統計モデルの予測の良さを表す量

# # 4.5.1 統計モデルの予測の良さ
# ## データ解析の目的
# + 観測される現象の背後のある「しくみ」の特定
# + 観測される現象の背後のある「しくみ」を近似的に代替しうる統計モデルの構築
# 
# ###  観測される現象の背後のある「しくみ」に近いか =「予測の良さ」
# ###  サンプルデータを説明できる = 「あてはまりの良さ」

# # 4.5.1 統計モデルの予測の良さ
# ## 予測の良さとは
# ### 同じデータ取得方法で異なるデータを得たとき、統計モデルはどれくらい正確に言い当てれるか
# ## 平均対数尤度とは
# ### 同じデータ取得方法で異なるデータセットを複数回あてはめた時の最大対数尤度の平均

# # 4.5.1 統計モデルの予測の良さ
# ## 例題 平均対数尤度の説明
# ### 1. パラメータ推定用データを作成する
# #### 真のモデル(β1=2.08=log(8))から50個の観測データを作る。
# 

# # 4.5.1 統計モデルの予測の良さ
# ## 例題 平均対数尤度の説明
# ### 2. 一定モデルのパラメータを推定する
# #### 50個の観測データで最尤推定し一定モデルのパラメータを推定する
# + 50個の観測データから推定された一定モデルのβ1が2.04とする
# ``` 
# logλi = β1 ・・・一定モデル
# 真のβ1=2.08
# 50個データを作ったときはβ1=2.04
# 最大対数尤度 LogL* = -120.6
# ```
# ![パラメータ推定用データを作成する](https://c4.staticflickr.com/6/5347/29289951403_5a5e518f0e.jpg "パラメータ推定用データを作成する")

# # 4.5.1 統計モデルの予測の良さ
# ## 例題 平均対数尤度の説明
# ### 3. 予測用データを作成する
# #### 真のモデル(β1=2.08=log(8))から50個の観測データから予測データセットを200個作成する
# ![真のモデル(β1=2.08)から50個データを作る。一定モデルのβ1が2.04となったとする ](https://c6.staticflickr.com/9/8053/29914071845_cafb83f916_b.jpg "真のモデル(β1=2.08)から50個データを作る。一定モデルのβ1が2.04となったとする")

# # 4.5.1 統計モデルの予測の良さ
# ## 例題 平均対数尤度の説明
# ### 4. 一定モデル(β1=2.04)で200個の予測データセットごとに最大対数尤度を求める
# ### 5. 200個の予測データセットごとの最大対数尤度から平均対数尤度を求める
# #### 200回最大対数尤度の平均が平均対数尤度
# ```
# 平均対数尤度 E(logL) = -122.9
# ```
# ![真のモデル(β1=2.08)から50個データを作る。一定モデルのβ1が2.04となったとする ](https://c6.staticflickr.com/9/8053/29914071845_cafb83f916_b.jpg "真のモデル(β1=2.08)から50個データを作る。一定モデルのβ1が2.04となったとする")

# # 4.5.2 最大対数尤度のバイアス補正
# ## 最大対数尤度と平均対数尤度の関係1
# 例題の最大対数尤度と平均対数尤度の関係
# + パラメータ推定用データの最大対数尤度 LogL'* = -120.06となる
# + 予測したときの　平均対数尤度 E(logL) = -122.9となる
# 
# ### 一定モデルの最大対数尤度が平均対数尤度 よりたかい
# #### 例題のパラメータはあてはまりの良さが過大評価されている
# ![パラメータ推定用データを作成する](https://c6.staticflickr.com/9/8120/29290473053_ac3c48aeb4.jpg "パラメータ推定用データを作成する")
# ※　パラメータ推定用データの最大対数尤度が平均対数尤度 を常に高くなるわけではない

# # 4.5.2 最大対数尤度のバイアス補正
# ## 最大対数尤度と平均対数尤度の関係2
# ```
# バイアス(b) = logL* -E(logL)
# 例題のバイアス(b) = -120.06 - -122.9 = 1.84
# ```
# ### １２回繰り返しを比較する
# ![最大対数尤度と平均対数尤度の計算を１２回](https://c7.staticflickr.com/9/8212/29623434590_7c0037d080.jpg "最大対数尤度と平均対数尤度の計算を１２回")

# # 4.5.2 最大対数尤度のバイアス補正
# ### 200回繰り返しを比較する
# ![最大対数尤度と平均対数尤度の計算を１２回](https://c3.staticflickr.com/9/8513/29623750090_03294eeb1b.jpg "最大対数尤度と平均対数尤度の計算を１２回")
# 例題のバイアスの標本平均は1.01

# # 4.5.2 最大対数尤度のバイアス補正
# ## 平均対数尤度E(logL)より、最大対数尤度logL*のほうが平均的に１ぐらい大きい

# # 4.5.2 最大対数尤度のバイアス補正
# ## 今回の例題の最大対数尤度と平均対数尤度の関係
# 最大対数尤度から平均対数尤度の推定値を求めるには、バイアスをひけばよい
# ```
# E(logL) = logL* -b
# bを1とすると
# E(logL) = logL* -1
# ```
# ## 数理統計学における最大対数尤度と平均対数尤度の関係
# ```
# E(logL) = logL* -k(パラメータ数)
# ```

# # 4.5.3 ネストしているGLM間のAIC比較
# ## 4.5.2でバイアス補正量(b)は1としたが、バラツキが多いのでは？
# ### ネストしている場合は、バラツキが小さくなる

# # 4.5.3 ネストしているGLM間のAIC比較
# #### ネストした関係のモデルの比較
# ```
# logλi = β1 ・・・一定モデル(k=1)
# logλi = β1 +  β2 * xi ・・・x モデル(k=2)
# ```
# ### 一定モデルとx モデルを比較する
# ####  ネストした関係一定モデルとx モデルと最大対数尤度の比較
# ```
# xモデルのlogL* - 一定モデルのlogL*  = 標本平均0.54　#xモデルのあてはまりが良くなる
# ```
# #### ネストした関係 一定モデルとx モデルと平均対数尤度の比較
# ```
# xモデルのE(logL) - 一定モデルのE(logL)  = 標本平均0.54　#xモデルの予測が悪くなる
# ```
# ![ネストしている一定モデルとxモデルの数値例](https://c5.staticflickr.com/9/8328/29837029092_b91375eb04.jpg "ネストしている一定モデルとxモデルの数値例")

# # 4.5.3 ネストしているGLM間のAIC比較
# ![ネストしている場合は、バラツキが小さくなる](https://c5.staticflickr.com/6/5697/29883862476_38ec805f34.jpg "ネストしている場合は、バラツキが小さくなる")
# ## 一定パラメータに一つ加えたモデルの概念的なまとめ
# + あてはまりのよさである最大対数尤度 log L*は0.5ぐらい増える
# + 予測能力である平均対数尤度E(log L)は0.5ぐらい減少する
# # 　　　　　　⇩
# ## logL* の平均バイアスは1増加する

# # 4.6 まとめ
# ## モデルを選択する方法
# ### 良いモデルとは
# #### 予測が良いものが良いモデルとかんがえるべき
# ## AICの使い方を理解する
# ## AICが正しい理由
# ### 追加したパラメータが最大対数尤度の統計的な変化量より大きいかどうかで選択できる

# # 参考資料
# 第４章を理解するために読んだ資料
#  - http://statistics.co.jp/reference/Toukeidatakaiseki_Nyumon/datakaiseki_nyumon10.pdf
#  http://www.it.mgmt.waseda.ac.jp/mi-tech/activity/data1/SelectionOfaStatisticalModel.pdf