担当 : 田中冬彦 2018 年 4 月 17 日 @ 統計モデリング 統計モデリング 第二回配布資料 文献 : A. J. Dobson and A. G. Barnett: An Introduction to Generalized Linear Models 3rd ed., CRC Press. 配布資料の PDF は以下からも DL できます. 短縮 URL http://tinyurl.com/lxb7kb8
第二回 ( 今回 ) 線形モデル 今後の予定 第三回一般化線形モデル Location 第四回ベイズ統計 ( 導入 ) Google map から転載 画像は You tube, 新唐人日本 2011 年 5 月 5 日付ニュース より * ニュースの内容の信ぴょう性には言及しません ( たとえば, 日本住血吸虫の場合, 経口感染でないことが実験で示されています ) ガンバ大阪ホームページより http://www2.gamba-osaka.net/stadium/ 第五回ベイズファクター Google map から転載
今日の内容 0. ( 統計の復習 ) 分布記号 & データの分類 1. 統計分析の流れ 2. 統計モデル 3. 線形モデル (4. グループ分けアンケート )
本日の主役 線形モデル ( 回帰モデル ) Y = α + β x + i i ε i i =1,2,, n ε i ~ N(0, 2 σ )
統計の復習 1 ~ 分布記号
分布記号 統計モデル = モデル式で表現 分布記号を使う 本講義でモデル式を使う理由 1. 統計 機械学習などのテキストで標準的に利用 2. WinBUGS, Stan などのツールで利用
分布記号の例 1 多項分布 M ( n; q1,, qk ) q 1 + q2 + + qk = 1 意味 : ツボの中に k 色の小さいボールを大量に入れる. その比率は q, q, 2, 1 q k n 個のボールを取り出す試行を考えるとき, 各色のボールの個数を X, X, 1 2, X k とする. これらは確率変数であり, 多項分布に従うことを以下のように記載. ( X, X 2,, X k ) ~ M ( n; q1,, q 1 k )
問 1: 練習してみよう! ツボに赤 (R) 青(B) 白(W) のボールを 5:3:2 の割合でいれてよく混ぜた. 100 個のボールを取り出す試行を考えるとき, 各色のボールの個数を X, X, X とする. R B W ( X, X, X R B W ) ~ M (100;0.5,0.3,0.2) 問 2: サイコロを 10 回ふって出た目の数を数える.(1~6 は 1/6 の確率で出る.) j の目が出る回数を X j ( X1, X 2, X 3, X 4, X 5, X 6) とする (j=1,2,3,4,5,6). ~ 1 M 10; 6, 1 6, 1 6, 1 6, 1 6, 1 6
分布記号の例 2 二項分布 意味 : Bin( n; q) 0 q 1 多項分布で k=2 ( 二色のボール ) を二項分布と呼ぶ. この場合, 片方の色のボールの個数のみに注目. ( 成功か失敗かの試行を n 回繰り返す ) X ~ Bin( n; q) 以下と同じ意味. ( X, Y ) ~ M ( n; q,1 q) X + Y = n
分布記号の例 3 正規分布 N( m, v) 意味 : 平均 m, 分散 v (>0) の正規分布 ( ガウス分布 ) X 確率変数が正規分布に従うことを以下のように記載. X ~ N( m, v) X, X, 2, X n 確率変数が同一の正規分布に独立に従うことを以下のように記載. (n 標本を独立に抽出, サンプリングする ) 1 i. i. d. X, X 2,, X ~ N( m, v) 1 n
練習してみよう! 問 3: 平均 162, 分散 25 の正規分布から 10 個の標本を抽出 X, X X 1 2,, 10 X, X,, X 1 2 10 i. i. d. ~ N(162,25) 補足 1. 分布記号のバリエーション X j j = 1,, n はすべて独立で以下の確率分布に従う X j ~ N( m, v) 2. 確率変数は通常 大文字だが 小文字で書いたり 混同して用いる
統計の復習 2 ~ データの分類
ここでの目標 データの形式と区別を理解 統計モデルとデータの形式は密接に関連 & モデルを紹介する際にデータのイメージと実例が思い浮かぶようにする
データの分類 (1/3) 変量 ( 変数 ) とは 1 変量データ x, x, 2, 1 x n 2 変量データ x, y ),( x, y ),,( x n, y ) ( 1 1 2 2 n k 変量データも同様に定義 (k 次元データとよぶことも ) n を標本数 ( サンプルサイズ ) とよぶ ( データサイズとよんだりすることもある ) 学食での摂取カロリー (kcal) 879 1047 760 779 845 ( 英語の点数, 統計の点数 ) (88, 90) (45, 78) (56, 100) (77, 85)
データの分類 (2/3) データの区分 ( 統計モデルを考える際の目安 ) 質的データ ( カテゴリカルデータ ) 量的データ ( 連続データ ) 名義尺度 順序尺度 間隔尺度 比率尺度 男 女 ( 性別 ) や職業など 〇 ( 評価 ) など ; 順序に意味があるが, 等間隔とは限らない 温度のように順序も間隔も意味があるが原点はどこでもよい 間隔尺度だが原点が定まっている. ( 重さ 長さなど ) * 参考 : 永田靖. 他著 : 多変量解析入門. サイエンス社, 1-1 節. 東京大学教養学部統計学教室編 : 統計学入門, 東京大学出版会, pp. 27-28.
モデリングする上での分類 データの分類 (3/3) 上限あり ある条件下での種子の発芽数 カウントデータ 上限なし その他のカテゴリカルデータ 交通事故件数 ; ツイート数 ; いいねの数 3 種類のラーメンの注文数 ( みそ しお とんこつ ) 正負をとる 温度 連続データ 正値のみ 製品の寿命
1. 統計分析の流れ
理想論 1. 分析課題 2. データ収集 3. データの統計分析 狭義にはここで 統計モデリング 4. 結論 実際には, 1,2,3,4 の順に進んで終了することはほとんどない!!
実際の所 例 1: まずはじめにデータありき IT 関係では大量のデータ 記録を保存 そこから 面白い関係を見つけ出してほしい ( むちゃぶりデータマイニング!) 例 2: 課題のすりかえ 分析したら 当初予定した結果が出なかった 1. 課題 も変更することに! 1,2,4 は完全に切り離して考えることはできない! 参考 : 松浦健太郎 : Stan と R でベイズ統計モデリング, 共立出版, Chap.3 統計モデリングを始める前に
2. 統計モデル
標本と母集団 統計学の復習 標本の例 : あるクラスの模試の点数 (72, 92, 91, 81, 73) 確率変数の実現値 ( 未知の分布 F から無作為に5つ取り出した値 ) と解釈 クラス B の受講者の点数分布 ( 仮想 72 92 91 81 73 0.00 0.01 0.02 0.03 0.04 0.05 F 0 20 40 60 80 100 この解釈により, 確率論と統計学が結びついた! 点数
標本と母集団 記法 X, X 2,, X i.i.d. ~ F 1 n 母集団 ( 分布 ) 観測される値の分布 ( 仮想的なものでもよい ) 統計の究極的な目標観測値から F が正確に把握できればよい 問題点 F の動く範囲は広すぎる ある程度, 分布の形を制限して考える ( だいたいの形がわかればよい )
統計モデルの設定 統計モデル : いくつかのパラメータで指定される確率分布の集合 pxθ ( ) θ 確率密度 / 確率関数確率分布の未知パラメータ x px ( θ)dx= 1, px ( θ) 0 px ( θ) = 1, px ( θ) 0 0.00 0.01 0.02 0.03 0.04 0.05 p( x θ ) クラス B の受講者の点数分布 ( 仮想 0 20 40 60 80 100 点数 記法 i.i.d. X,, ~ ( ) 1 Xn pxθ
最初の統計モデリング シチュエーション 2 種類の方法 A, B で金を回収 [g]. 廃棄携帯の基盤ひと山あたりの回収量が A, B で以下のようになった. 基本的な統計量 A: 73, 72, 66, 80, 75 B: 71, 67, 68, 57, 68, 75, 60, 69 全体の平均 69.3 全体の分散 38.4 Aの平均 73.2 Aの分散 25.7 Bの平均 66.9 Bの分散 33.5 なんとなく A の方が回収量が多い?( 金なので 差は無視できない )
最初の統計モデリング 統計モデルの設定二つとも連続値 とりあえず, 正規分布からの標本と仮定分散は等しい ( 解析を簡単化する仮定 ) モデル式 ( 独立な2 変量ガウスモデル ) i.i.d. X, X,, X ~ N(, v) 1 2 5 i.i.d. Y, Y,, Y ~ 1 2 8 N( µ A µ B, v) µ,σ 注意 : 統計モデルのパラメータは, p, q, f, t, など何を用いてもよい. ただし, 異なるものはのように区別すること. µ A, µ B
* 計算公式は省略 ( 統計のテキストに掲載 ) 最初の統計モデリング モデルパラメータ ( 母数 ) の推定値 * ˆ = 73.2, µ A ˆ µ B = 66.9 vˆ = 30. 7 Ambition of TKK 注 : パラメータの推定量 ( 値 ) はハットをつける 可視化の例 パラメータの推定値を代入して分布を眺める Population 0.00 0.02 0.04 0.06 0.08 40 50 60 70 80 90 100 yields A の方が B より回収量が多め ( 本来はこの後, t 検定 )
ここまでのまとめ 統計モデリングの基本的な考え方 1. データ ( 数値 ) の背後に母集団分布を想像 2. 母集団分布を統計モデルで表現 パラメータ推定 ( 点推定 ) や信頼区間 仮説検定 予測 課題が先か手法が先か仮説検定 ( 統計手法 ) を知っていると それに応じた課題設定が可能 分析課題 : 方法 A, B で回収量に差があるか? 分析手法 : ( ガウスモデルでの ) 平均の差の仮説検定
練習してみよう! O 大学 ( 数千人規模 ) から無作為に 100 人の学生を選び出し, A, B,C 三択のアンケートを行い以下のような結果を得た. 分析のためのモデル式を設定しなさい. A: 85 B: 13 C: 2 合計 : 100 1) O 大学で A, B,C と回答する人たちの割合を以下のように表す ( 分析者にとって未知のパラメータ ) q, q, q, q + q + q = 1 ( ) A B C A B C
さらに 無作為に選んだ学生がアンケートでA, B, Cと回答する人数をそれぞれ X A, XB, XC とすると これは確率変数と考えることができる 2) 無作為に一人を選んだ場合, その人が B と回答する確率はいくらか? パラメータで書きなさい. q = = = = B ( 0, 1, 0) P X X X A B C 3) 無作為に三人を選んだ場合, 一人が A, 二人が B と回答する確率はいくらか? パラメータで書きなさい. ( 1, 2, 0) P X = X = X = = A B C 3qq A B 2
4) 同様に 100 人選んだ場合, それぞれの回答が x, x, x, x + x + x = 100 ( ) A B C A B C のようになる確率はいくらか? P X = x, X = x, X = x = ( ) A A B B C C 100! x x q q q x! x! x! A B C A B x C A B C 5) 4) の結果を分布記号を使ってモデル式で表しなさい A, ~ M(100; q, q, q ) B C ( X X, X ) A B C 参考 学部 1 年で習う母比率の差の仮説検定の公式は上のような 3 項モデルで導出している.
モデルを設定する理由 補足 1. 母集団分布の正確な形状は知り得ない, だいたいの形で十分 2. 実験結果から分布の形状が既知の場合, 正当化できる (*) 3. 仮説検定や信頼区間 初等統計では, このことをあまり表に出さない ベイズ分析で必要 * 精密科学 / 実験科学の状況だが, 本講義ではあまり考えないシチュエーション, 思想的な注意点 1. 通常, モデルに正解はない / 検証のしようがない物理などの自然科学 2. 独立同一性 (i.i.d.) の仮定も含め 作業仮説 3. よいモデル は目的 課題依存
3. 線形モデル
ここでの目標 ある変数を別の変数で説明するモデルを提案 & モデルパラメータの推定 注 : ここでの例は分析課題は提示しない
回帰分析 (B-2/C-2 資料より ) 例題 : みずほの部屋探し O 大学新入生のみずほさんは賃貸情報をネットで検索. 以下のようなデータを得ました. 豊中キャンパス近くの賃貸物件 (1K) 最寄り駅からの距離 ( 徒歩 ): 3 5 6 10 11 17 一カ月の賃料 ( 万円 ): 8 7.3 6.2 4 4.2 3.5 傾向をみるため, 横軸に距離, 縦軸に賃料をとりプロット ( 点を打つ )
データのプロット ワンポイント ペアになっている 2 変量データは プロットしておおまかな傾向をつかむ! R プログラム例 x <- c(3, 5,6, 10, 11, 17); y <- c(8, 7.3, 6.2, 4, 4.2, 3.5); plot(x,y, pch=18, col=2, xlim=c(0, 20), ylim=c(0, 10), main="kaiki", xlab="min Walk", ylab="10^4 YEN"); abline(h=0, lty=2, col="gray"); # hori line abline(v=0, lty=2, col="gray"); # vert line 10^4 YEN 0 2 4 6 8 10 Kaiki 0 5 10 15 20 Min Walk なんとなく右肩下がりの傾向が見える
説明変数と目的変数 説明変数と目的変数 ( データのばらつきはいったん無視 ) 簡単な関数 f で変数に以下のような関係が期待される時 y x 説明変数 y 目的変数 f (x) y 10^4 YEN 0 2 4 6 8 10 Kaiki とよぶ. 今回はx, y は1 次元 (1 変量 ) のみ扱う. ( 因果関係が既知の ) 統計モデリング y f ( x 1,, x k ) この f をうまく与える ( モデル化 ) のがひとつの目標 0 5 10 15 20 Min Walk x
統計モデルの導入 統計モデルの設定 y f (x)?? なんとなく右肩下がり とりあえず, f として直線 ( 一次式 ) を仮定 f ( x) = α + β x ε i : = yi ( α + β xi ) とりあえず, 平均 0 の正規分布を仮定 分散は等しい ( 解析を簡単化する仮定 ) 本来のモデル式 ε,, i.i.d. ~ N(0, σ 2 ), 2 6 ε ε 1
線形モデル ( 回帰モデル ) 線形モデル 通常は, 以下のような形で記載 ( f(x) の形を明示 ) Y = α + β x + i i ε i i =1,2,,6 2 ε ~ N(0, σ ) i x: 最寄駅からの距離 ( 分 : 徒歩換算 ), y: 一か月の家賃 ( 万円 ) モデルのパラメータ, α, β; σ パラメータは最尤推定法などで推定 (R コマンドでできる ) ˆ α = 8.5 ˆ β = 0. 34 推定値を代入した f(x) ( 回帰直線という ) y = ˆ α + ˆ β x = 8.5 0. 34x 2
回帰直線 R プログラム例 ( 回帰分析 ) x <- c(3, 5,6, 10, 11, 17); y <- c(8, 7.3, 6.2, 4, 4.2, 3.5); res <- lm(y~x); ahat <- res$coefficients[1]; bhat <- res$coefficients[2]; R プログラム例 ( 回帰直線 ) plot(x,y, pch=18, col=2, xlim=c(0, 20), ylim=c(0, 10), main="kaiki", xlab="min Walk", ylab="10^4 YEN"); abline(h=0, lty=2, col="gray"); # hori line abline(v=0, lty=2, col="gray"); # vert line abline(a=ahat, b=bhat); 10^4 YEN 0 2 4 6 8 10 Kaiki y = ˆ α + ˆ β x = 8.5 0. 34x 0 5 10 15 20 Min Walk
ここまでのまとめ 今の例について 1K の家賃は 最寄駅からの距離 ( 徒歩換算 ) が増えるほど 減少する傾向がみてとれた だいたい一次式に従っている より踏み込んだ分析に向けて あてはまりのよさも議論 ( 仮説検定 ) 最寄駅からの距離で だいたいの家賃を予測
さまざまな拡張のアイディア 1. 目的変数 Y の期待値を x で表現 ( 確率の要素なし!) µ : = EY [ ] = f( x) i i i 注 : 確率変数の期待値 = α + β x i たとえば, 多項式回帰など. f ( x) = α + β x + γ x 2 2.Y の分布のモデル Y N x i 2 ~ ( α + β i, σ ) 参考 : 解釈無視で, x, y を log, べき乗で変換する方法もある
複数のモデルがある場合 * 本講義では, AIC, BIC などの情報量規準を機械的に使用してよい 赤池情報量規準 (Akaike Information Criterion) を用いたモデル選択 AIC = 2 L( θ ˆ) + 2 p L (θ ) 尤度関数 θˆ 最尤推定での推定値 p パラメータ数 尤度関数の最大値 L ( ˆ) θ = max L( θ ) ( 最尤推定は尤度関数を最大化 ) とパラ θ メータ数からAICを計算 ; 相対的にAICの小さいモデルを選ぶ グループタスクでは統計研究室の学生に聞こう!