スライド 1

2017年4月18日＠統計モデリング

担当：田中冬彦

統計モデリング

第二回 配布資料

文献

:

A. J. Dobson and A. G. Barnett:

An Introduction to Generalized Linear Models 3rd ed.,

CRC Press.

配布資料の

PDFは以下からもDLできます.

短縮URL http://tinyurl.com/lxb7kb8

今後の予定

Google map から転載 Location

第二回（今回） 線形モデル

第三回 一般化線形モデル

*ニュースの内容の信ぴょう性には言及しません (たとえば, 日本住血吸虫の場合, 経口感染でないことが実験で示されています。) 画像はYou tube, 【新唐人日本2011年5月5日付ニュース】より

第四回 ベイズ統計（導入）

ガンバ大阪ホームページより http://www2.gamba-osaka.net/stadium/ Google map から転載

第五回 ベイズファクター

今日の内容

０． （統計の復習） 分布記号＆データの分類

１．統計分析の流れ

２．統計モデル

３．単回帰の統計モデル

４．線形モデル

（５．グループ分けアンケート）

本日の主役

線形モデル（単回帰モデル）

6

,

,

2

,

1 

=

i

)

,

0

(

~

σ

2

ε

_i

N

i

i

i

x

Y

=

α

+

β

+

ε

統計の復習１

～分布記号

分布記号

統計モデル ＝ モデル式で表現

↑ 分布記号を使う

本講義でモデル式を使う理由

１． 統計・機械学習などのテキストで標準的に利用 ２． WinBUGS, Stan などのツールで利用

分布記号の例１

ツボの中に

k

色の小さいボールを大量に入れる. その比率は 多項分布 意味： k

q

q

q

₁

,

₂

,



,

1

2 1

+

q

+

+

q

k

=

q



n

個のボールを取り出す試行を考えるとき, 各色のボールの個数を k

X

X

X

₁

,

₂

,



,

とする. これらは確率変数であり, 多項分布に従うことを以下のように記載.

)

,

,

;

(

~

)

,

,

,

(

X

₁

X

₂



X

_k

M

n

q

₁



q

_k

)

,

,

;

(

n

q

₁

q

_k

M



ツボに赤（R）・青(B)・白(W)のボールを、５：３：２ の割合でいれてよく 混ぜた. １００個のボールを取り出す試行を考えるとき, 各色のボール の個数を

練習してみよう！

問１：

~

)

,

,

(

X

_R

X

_B

X

_W

W B R X X X , , とする.

)

2

.

0

,

3

.

0

,

5

.

0

;

100

(

M

問２： サイコロを１０回ふって出た目の数を数える.（１～６は1/6 の確率で出る.） j の目が出る回数を とする (j=1,2,3,4,5,6). X j

~

)

,

,

,

,

,

(

X

₁

X

₂

X

₃

X

₄

X

₅

X

₆













6

1

,

6

1

,

6

1

,

6

1

,

6

1

,

6

1

;

10

M

分布記号の例２

二項分布 意味： 多項分布で k=2 (二色のボール)を二項分布と呼ぶ. この場合, 片方の色のボールの個数のみに注目. (成功か失敗かの試行を n 回繰り返す)

1

0

≤ q

≤

以下と同じ意味.

)

1

,

;

(

~

)

,

(

X

Y

M

n

q

−

q

)

;

(

n

q

Bin

)

;

(

~

Bin

n

q

X

n

Y

X

+

=

分布記号の例３

正規分布 意味：

)

,

(

N

~

,

,

,

₂

1

X

X

m

v

X



_n

確率変数 が正規分布に従うことを以下のように記載.

)

,

(

N

m

v

平均 m, 分散 v (>0) の正規分布（ガウス 分布）

X

)

,

(

N

~

m

v

X

確率変数 が同一の正規分布に独立に従うこと を以下のように記載. （n 標本を独立に抽出, サンプリングする） n

X

X

X

₁

,

₂

,



,

. . . di i

平均 162, 分散 25 の正規分布から１０個の標本 を抽出

練習してみよう！

問３： 10 2 1, X , , X X 

~

,

,

,

₂

₁₀

1

X

X

X



i. di. .

N

(

162

,

25

)

１．分布記号のバリエーション

)

,

(

N

~

m

v

X

_j n j =1 , , はすべて独立で以下の確率分布に従う j

X

補足

２．確率変数は通常、大文字だが、小文字で書いたり、混同して用いる

統計の復習２

～データ分類

ここでの目標

１．世の中（統計の本）には

色々な形式のデータ

がある

ことを理解

２．「用語」を暗記する必要なし！

→ モデルを紹介する際に

データのイメージと実例

が

思い浮かぶようにする

変量（変数）とは

ｋ変量データも同様に定義 (k次元データとよぶことも) (英語の点数, 統計の点数) 88, 90 45, 78 56, 100 １変量データ 体脂肪率の減少量 -0.08 10.47 10.87 -12.28 n x x x₁, ₂,, ２変量データ (x₁, y₁),(x₂, y₂),,(x_n, y_n) n を標本数 (サンプルサイズ)とよぶ （データサイズとよんだりすることもある）

データの分類(1/3)

データの分類(2/3)

データの区分

量的データ （連続データ） 質的データ （カテゴリカルデータ) ・名義尺度 ・間隔尺度 男、女（性別）や職業など ◎、〇、△、×（評価）など；順 序に意味があるが, 等間隔とは 限らない *参考： 永田 靖. 他 著: 多変量解析入門. サイエンス社, 1-1節. 東京大学教養学部統計学教室編: 統計学入門, 東京大学出版会, pp. 27-28. 温度のように順序も間隔も意味 があるが原点はどこでもよい ・比率尺度 ・順序尺度 間隔尺度だが原点が定まって いる. (重さ、長さなど)

モデリングする上での分類

連続データ カウントデータ ・上限あり ・正負をとる ある条件下での種子の発芽数 交通事故件数 ３種類のメニューの注文数 (みそ、しお、とんこつ) ・正値のみ ・上限なし 温度 その他のカテゴリカルデータ 製品の寿命

データの分類(3/3)

理想論

実際には, 1,2,3,4 の順に進んで終了することはほとんどない！！

← 狭義にはここで「統計

モデリング」

１．分析課題

３．データの統計分析

４．結論

２．データ収集

実際の所

IT関係では大量のデータ・記録を保存 → そこから、面白い関係を見つけ出してほしい （むちゃぶりデータマイニング！）

例１： まずはじめにデータありき

例２： 課題のすりかえ

分析したら、当初予定した結果が出なかった → 「１．課題」も変更することに！ 実際には, １，２，４は完全に切り離して考えることはできない！ 参考： 松浦健太郎: StanとRでベイズ統計モデリング, 共立出版, Chap.3 「統計モデリングを始める前に」

この解釈により, 確率論と統計学が結びついた！

標本の例：

あるクラスの模試の点数 (72, 92, 91, 81, 73)

確率変数の実現値

(未知の分布 F から無作為に５つ取り出した値)

と

_解釈

0 20 40 60 80 100 0. 00 0. 01 0. 02 0. 03 0. 04 0. 05 クラス B の受講者の点数分布(仮想 点数

F

72 92 91 81 73

標本と母集団

~

,

,

,

₂ 1

X

X

n

X



i.i.d.

F

標本と母集団

記法

母集団(分布)

観測される値の分布

問題点

Fの動く範囲は広すぎる

→ ある程度, 分布の形を制限して考える

0 20 40 60 80 100 0. 00 0. 01 0. 02 0. 03 0. 04 0. 05 クラス B の受講者の点数分布(仮想 点数

)

|

(

x

θ

p

統計モデル

_{： いくつかのパラメータ で指定される}

確率分布の集合

統計モデルの設定

ただし, メジャーな分布記号を用いることも多い 確率分布の未知パラメータ

θ

i.i.d.

)

|

(

~

,

,

1

y

p

y

θ

y



_n

記法(一例)

)

|

(

y

θ

p

_∫

_p

₍

_y

_|

_θ

₎

_d

_y

₌

₁

_,

_p

₍

_y

_|

_θ

₎

_≥

₀

0

)

|

(

,

1

)

|

(

=

≥

∑

p

y

θ

p

y

θ

最初の統計モデリング

２種類の方法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の方が回収量が多い？

（金なので、差は無視できない）

最初の統計モデリング

二つとも連続値

→ とりあえず,

_{正規分布からの標本と仮定}

モデル式 （独立な２変量ガウスモデル）

~

,

,

,

₂

₅

1

X

X

X



i.i.d.

)

,

(

N

µ

_A

v

8

2

1

,

Y

,

,

Y

Y



_~

i.i.d.

)

,

(

N

µ

_B

v

統計モデルの設定

σ µ, 注意：統計モデルのパラメータは, p, q, f, t, など何を用いてもよい. ただし, 異なるものは のように区別すること. µA,µB

分散は等しい（解析を簡単化する

_仮定)

最初の統計モデリング

モデルパラメータ(母数)の推定値*

可視化の例

,

2

.

73

ˆ

_A

=

µ

*計算公式は省略（統計のテキストに掲載）

パラメータの推定値を代入して

分布を眺める

9

.

66

ˆ

_B

=

µ

v

ˆ

=

30

.

7

Aの方がBより回収量が多め（本来はこの後, t 検定）

40 50 60 70 80 90 100 0. 00 0. 02 0. 04 0. 06 0. 08 Ambition of TKK yields P opul at ion 注：パラメータの推定量(値)はハットをつける

ここまでのまとめ

１． データ（数値）の背後に母集団分布を想像

２． 母集団分布を統計モデルで表現

→ パラメータ推定（点推定）や信頼区間、仮説検定、予測

統計モデリングの基本的な考え方

課題が先か手法が先か

分析手法： （ガウスモデルでの）平均の差の仮説検定 分析課題： 方法A, Bで回収量に差があるか？

仮説検定（統計手法）を知っていると、それに応じた課題

設定が可能

練習してみよう！

O大学（数千人規模）から無作為に１００人の学生を選び出し,

A, B,C三択のアンケートを行った.

三項分布の記号を使ってモデル式を書きなさい.

A: 85

B: 13

C: 2

合計: 100

モデル式

(

X

_A

,

X

_B

,

X

_C

)

~

M

(

100

;

p

,

q

,

r

)

補足

思想的な注意点

１． たいていの場合, 正解はない/ 検証のしようがない

２． 独立同一性(i.i.d.)の仮定も含め、作業仮説

３． 「よいモデル」は目的・課題依存

1. 母集団分布の正確な形状は知り得ない,

形状に興味はない（誤差モデル

→ 分散が重要）

2. 実験結果から分布の形状が既知の場合, 正当化できる(*)

3. 仮説検定や信頼区間、ベイズ分析で必要

モデルを設定する理由

*精密科学/実験科学の状況だが, 本講義ではあまり考えないシチュエーション,

ここでの目標

ある変数を別の変数で説明するモデルを提案

＆ モデルパラメータの推定

回帰分析（B-2/C-2資料より）

O大学 新入生のみずほさんは賃貸情報を ネットで検索. 以下のようなデータを得ました.

例題： みずほの部屋探し

→ 傾向をみるため, 横軸に距離, 縦軸に賃料をとりプロット(点を打つ)

最寄り駅からの距離 (徒歩)： 3 5 6 10 11 17

一カ月の賃料 (万円)：

8 7.3 6.2 4 4.2 3.5

豊中キャンパス近くの賃貸物件（１K）

データのプロット

0 5 10 15 20 0 2 4 6 8 10 Kaiki Min Walk 10^ 4 Y E N 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 R プログラム例

ペアになっている２変量データは

プロットしておおまかな傾向をつかむ！

ワンポイント

説明変数と目的変数

0 5 10 15 20 0 2 4 6 8 10 Kaiki Min Walk 10^ 4 Y E N （データのばらつきはいったん無視） 簡単な関数 f で変数に以下のような関係が期待される時

説明変数と目的変数

)

( x

f

y

≈

y

目的変数

x

説明変数 とよぶ.

(因果関係が既知の) 統計モデリング

講義では目的変数は１次元(1変量)のみ扱う.

)

,

,

(

x

₁

x

_k

f

y

≈



この f をうまく与える（モデル化）のがひとつの目標

x

y

統計モデルの導入

)

(

:

_{i i} i

y

α

β

x

ε

=

−

+

統計モデルの設定

??

)

( x

f

y

≈

なんとなく右肩下がり

→ とりあえず,

f として直線（一次式）を仮定

モデル式

~

,

,

,

₂

₆

1

ε

ε

ε



i.i.d.

_N

₍

₀

_,

σ

2

₎

分散は等しい（解析を簡単化する

_仮定)

→ とりあえず, 平均０の正規分布を

仮定

x

x

f

(

)

=

α

+

β

6 , , 2 , 1  = i x: 最寄駅からの距離(分：徒歩換算), y:一か月の家賃 (万円)

線形モデル

)

,

0

(

~

σ

2

ε

_i

N

i

i

i

x

Y

=

α

+

β

+

ε

モデルのパラメータ,

α

,

β

;

σ

2

パラメータは最尤推定法などで推定 (Rコマンドでできる)

x

x

y

=

α

ˆ

+

β

ˆ

=

8

.

5

−

0

.

34

線形モデル（回帰モデル）

通常は, 以下のような形で記載 ( f(x) の形を明示)

5

.

8

ˆ

=

α

β

ˆ

=

−

0

.

34

推定値を代入した f(x) (回帰直線という)

回帰直線

0 5 10 15 20 0 2 4 6 8 10 Kaiki Min Walk 10^ 4 Y E N 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);

R プログラム例 （回帰直線）

x x

ここまでのまとめと補足

・１Kの家賃は、最寄駅からの距離（徒歩換算）が増えるほ

ど、減少する傾向がみてとれた。

・だいたい一次式に従っている

今の例について

より踏み込んだ分析に向けて

・あてはまりのよさも議論（仮説検定）

・最寄駅からの距離で、だいたいの家賃を予測

1次式でうまくいかない場合

・解釈無視で, x, y をlog, べき乗で変換

・多項式回帰など.

2

)

(

x

x

x

f

=

α

+

β

+

γ

ここでの目標

データをあれこれ分析してから

逆に課題を設定する流れを理解

Birthweight vs Gestational Age

(余計なものは取り除いてある) 胎内にいた期間 [週], 出生時の体重 [g], 男児(b)/女児(g) 男児, 女児 ともに標本サイズは １２ずつ

データ例

1 40 2968 b 2 38 2795 b 3 40 3163 b 4 35 2925 b 5 36 2625 b 6 37 2847 b 7 41 3292 b 8 40 3473 b 9 37 2628 b 10 38 3176 b 11 40 3421 b 12 38 2975 b 13 40 3317 g 14 36 2729 g 15 40 2935 g 16 38 2754 g 17 42 3210 g 18 39 2817 g 19 40 3126 g 20 37 2539 g 21 36 2412 g 22 38 2991 g 23 39 2875 g 24 40 3231 g

35 36 37 38 39 40 41 42 2400 2600 2800 3000 3200 3400 Chap. 2 Age W ei ght b g 定量的な確認： たとえば相関係数の計算

ρ

=

0

.

744

データのプロット

見てわかること

24 , , 2 , 1  = i x: 胎内にいた期間(週), y:出生時の体重 i i i

x

Y

E

[

]

=

µ

=

α

+

β

)

,

0

(

~

σ

2

ε

_i

N

統計モデルの導入

i

i

i

x

Y

=

α

+

β

+

ε

モデルのパラメータ,

α

,

β

;

σ

2

パラメータを最尤推定法によって推定し以下の直線を引いてみる

x

Y

=

α

ˆ

+

β

ˆ

線形モデル （まずは, 男女の区別なし, ２４のデータと考えて分析）

回帰直線を引いてみる

データ > dat2

AGE WEI TYPE 1 40 2968 b 2 38 2795 b ・・・・ 23 39 2875 g 24 40 3231 g 線形回帰

> dat2.res <- lm(WEI~ AGE, data= dat2);

35 36 37 38 39 40 41 42 2400 2600 2800 3000 3200 3400 Chap. 2 Age W ei ght b g 回帰直線

x

y

=

α

ˆ

+

β

ˆ

α

ˆ

=

−

1484

,

β

ˆ

=

115

.

5

R プログラム例

データを男児(j=1), 女児(j=2)に分けて分析 線形モデル

)

,

0

(

~

σ

2

ε

_ji

N

ji ji j j ji

x

Y

=

α

+

β

+

ε

回帰直線

x

y

=

α

ˆ

_j

+

β

ˆ

_j

112

ˆ

,

1269

ˆ

₁

=

−

β

₁

=

α

35 36 37 38 39 40 41 42 2400 2600 2800 3000 3200 3400 Chap. 2 Age W ei ght b g

130

ˆ

,

2142

ˆ

₂

=

−

β

₂

=

α

線形回帰はよさそうだが, 男児と女児で分けて考えた方がいいのか？ 分析の課題

男児・女児に分けて推定すると

２種類の線形モデルで仮説検定（一般的な形） H1: 傾きが異なる線形モデル

)

,

0

(

~

σ

2

ε

_jk

N

jk jk j j jk

x

Y

=

α

(1)

+

β

+

ε

仮説検定によるモデル選択 (1/2)

H０: 傾きは等しい線形モデル jk jk j jk

x

Y

=

α

(0)

+

β

+

ε

ε

_~

₍

₀

_,

σ

2

₎

N

jk 検定統計量：

∑

−

−

=

k j jk j jk

x

Y

S

, 2 ) 0 ( 0

:

α

ˆ

β

ˆ

=

∑

− − k j jk j j jk x Y S , 2 ) 1 ( 1: αˆ βˆ ) 2 ( , 1 1 1 0

_~

1

)

2

(

− −

−

−

−

K J J

F

J

K

J

S

S

S

12 ; 2 = = K J K k J j =1,..., ; =1,...,

２種類の線形モデルで仮説検定 H1: 傾きが異なる線形モデル H０: 傾きは等しい線形モデル VS (有意水準 0.05)仮説H0が正しいとすると, 検定統計量の値は0から4.35程度 におさまるはず → 実際に（がんばって）計算すると・・・

35

.

4

19

.

0

1

)

2

(

1 1 0

=

<<

−

−

−

J

K

J

S

S

S

（

H0 は棄却されないため）傾きが異なるとは言えない

仮説検定によるモデル選択 (2/2)

ここまでのまとめと補足

・男女ともに, 出生時の体重は胎内にいる期間に対し,

一次式にしたがって増える

・(切片は違うが) 男女間で傾きに有意な差は認められない

今の例について

複数のモデルがある場合

（仮説検定で考えるケースは稀）

本講義では, AIC, BICなどの情報量規準を機械的に使用

してよい

p

L

AIC

=

−

2

(

θ

ˆ

)

+

2

赤池情報量規準

(Akaike Information Criterion)

を用いたモデル選択

尤度関数の最大値 (最尤推定は尤度関数を最大化)とパラ メータ数からAICを計算; 相対的にAICの小さいモデルを選ぶ p パラメータ数 ) ( max ) ˆ (θ θ θ L L = ) (θ L 尤度関数 θˆ 最尤推定での推定値