多変量解析とは

Excel入門

ロジスティック回帰

教科書のデータ

表7.2 かぶと虫の死亡データ ガス濃度はlog10CS2mgl-1

ガス濃度 かぶと虫の数 死亡数 x n y 1.6907 59 6 1.7242 60 13 1.7552 62 18 1.7842 56 28 1.8113 63 52 1.8369 59 53 1.861 62 61 1.8839 60 60 出典

ロジスティック回帰

表7.2から、ガス濃度毎の死亡割合を計算して図にしました。

このデータに基づいて、死亡確率をガス濃度の関数として表したいと思います。但し確率は0から1の間ですので、直線や多項式を用いると、ガス濃度によっては0から1の範囲をはみ出てしまうので不適です。そこで0から1の範囲をはみ出さないモデルとして次の3つがしばしば使われます。
表の上からi番目の行は、ガス濃度がx_iであるときに、かぶと虫n_i匹の中のx_i匹が死んだということを表しますので、その濃度で各々のかぶと虫が死ぬ確率をp_iとすると、死亡数の確率分布は二項分布B(n_i,p_i)になります。このp_iをx_iのどのような関数で表すか、という問題を考えます。二項分布の確率関数p(x)=_nC_xp^x(1-p)^n-xで、Excelでは_nC_x=COMBIN(n,x)です。実はこの項にはp(x)が含まれないので省いも構いません。

モデル1:ロジスティックモデル

死亡確率はp(x)=exp(α+βx_i)/(1+exp(α+βx_i)))である。

モデル2:プロビットモデル

死亡確率はp(x)=Φ((x_i-μ)/σ)である。但しΦは標準正規分布の分布関数で、Excelでは Φ((x-μ)/σ))=NORM.DIST(x,μ,σ,1)で書けます。

モデル3:極値モデル

死亡確率はp(x)=1-exp(-exp(α+βx_i))である。

まずロジスティックモデルで確率を計算しましょう。α、βは最初は分からないので、xの係数βは0にして、αはp(x)=1/2になるようにα=0にします。

次に対数尤度の計算です。二項分布B(n,p)の確率関数は_nC_xp^x(n-p)^n-xでした。対数尤度はlog_nC_x+x_ilog p_i+(n_i-x_i)log (1-p_i)をiを1から8まで全部足したものになります。

F12セルの対数尤度を最大にするようなF1、F2セルの値を、ソルバーを用いて求めてください。

モデルによって推定した死亡確率と、実際の死亡割合を比較するとこのようになります。

点は実際の死亡割合、曲線がモデルです。

プロビットモデルに関しても同様にまず確率を求めて、それから対数尤度を計算してみましょう。μは最初はデータの中ごろということで1.8、σはとりあえず1にしてソルバーを使ってください。

すると、分母のσが小さくなりすぎてエラーになりました。

エラーにならないように変数の範囲に制約をつけましょう。

σ=0.02…でエラーになっているので、σ≧0.03にしましょう。制約条件(U)の追加(A)をクリックして、制約式を入力します。

もう一度ソルバーを実行してみます。

今度はエラーになりませんでしたし、σ=0.05…ですので制約の端でもありません。もしもソルバーの解が制約の端である0.03と一致していたら、ソルバーの結果が制約式に依存してしまっているので、制約式を少し広げてσ≧0.025などを試すべきですが、今回は端ではないのでこれで大丈夫です。

最後に極値モデルに関しても最尤推定量を求めてください。α、βの初期値はロジスティックモデル同様α=β=0にします。今度もソルバーでエラーが出ました。

今回は死亡確率が1になっているのは最後の行だけで、ここは60匹全員が死亡しているので尤度としては1⁶⁰で良いのですが、公式にそのまま代入すると_nC_xp^x(1-p)^n-xの中の(1-p)^n-xの部分が0⁰となり、数学的には1なのですがExcelは底が0だとエラーになってしまいます。エクセルがエラーを起こす(1-p)^n-xの部分は、最後の行ではどうせ1を掛けているだけなので削除してしまっても計算式は変わりません。

今回はソルバーで計算することが出来ました。

レポート課題

それぞれのモデルの最尤推定量を求めて、どのモデルが一番良いか選択したエクセルファイルを提出してください。

コンペティション

夏休み中に参加できるコンペティションの一覧です。

• 第6回 データビジネス創造コンテスト
• 第7回スポーツデータ解析コンペティション
• マーケティング分析コンテスト2017
• 電力需要予測コンテスト
• 地方創生政策アイデアコンテスト2017

品質管理

品質管理においてベイズの定理を使って計算しました。講義ではこの他にも別の先生の資料も使いました。

ブートストラップ法

前回までの解答をパワーポイントで説明した後、板書でブートストラップについて説明しました。

ポアソン分布

5月22日の講義に引き続いて、Excelのソルバーを用いて最尤推定量を求めます。5月22日に欠席してソルバーの準備をしていない人は5月22日の講義のページを見て準備してください。
注:5月22日の講義の時と別のパソコンを使っている人は、改めてソルバーの準備が必要かも知れません。

教科書のデータ

これまでの講義で何度か使った、都市部と農村部の慢性病の数のデータです。 一般化線形モデル入門22ページ表2.1のデータです。

都市群 0 1 1 0 2 3 0 1 1 1 1 2 0 1 3 0 1 2 1 3 3 4 1 3 2 0

農村群 2 0 3 0 0 1 1 1 1 0 0 2 2 0 1 2 0 0 1 1 1 0 2

出典

これらのデータが、どんな母集団分布に従って観測されたのか、特に都市群と農村群で同じ分布と考えるのが妥当なのか、違うと考えるのが妥当なのか、考えます。
母集団分布を寸分の狂いも無く推定するには無限個のデータが必要ですのでそれは諦めて、ポアソン分布を用いて推定してみます。他にももっと適当な分布があるかも知れません。
ポアソン分布の確率関数は f(x|λ)=λ^xe^-λ/x! です。そしてそのグラフをいくつかのλに関して書いてみると

になります。今日はまず都市群のデータ

に一番良く当てはまるλの値を最尤法を用いて求めてみましょう。

パラメータのベクトルをθとすると、対数尤度関数は∑_α=1ⁿlog(f(x_α|λ))です。ポアソン分布の場合はパラメータはλ一つですから、θ=λです。ポアソン分布の場合は対数尤度関数をλで微分して=0とおいた式を解くことでλをx₁x₂,…,x_nを用いて表すことも出来ますが、今回は解けない場合に計算機を使ってどうやって求めるか、という問題を考えていますので、Excelで尤度関数の値を求めて、それを最大にするλを探しましょう。

まず都心部のデータをExcelに貼り付けます。説明を合わせるためにA1セルから貼り付けてください。データはA2からA27に入ります。次にλの値を入力します。後で変化させますがまずはλ=1にしましょう。B1セルにλ、C1セルに1と書いてください。
そして対数尤度関数はlog(f(x_α|λ))を足したものですので、B2セルに、この式のxのところにA2セルの値が、λのところにC1セルの値が来るように式を書きます。但し、この式を後でB3以降にコピーしますので、C1はC$1と書いて、下へコピーしても番号が変わらないようにしてください。
ポアソン分布の確率関数の式はEXCELでは =POISSON.DIST(x, λ, 0) です。EXCEL2007以前では =POISSON(x, λ, 0) と書きます。また、底がeである自然対数はLN( )です。

練習

D1セルに都市群のデータの対数尤度関数を計算し、C1セルの値を変化させて、対数尤度が最大になるλの値をソルバーを用いて探してください。
また、対数尤度関数をλで微分して=0とおいた式を解いて、λの最尤推定量をx₁x₂,…,x_nを用いた式で書き、その式に都市群のデータを代入した値と比較してください。

次に二通りの統計モデルを考えて、最尤推定量を求めます。
都市群と農村群で共通のパラメータを用いるモデルと、別のパラメータを用いるモデルです。

まず農村部のデータを貼り付けます。ひとまとめにして扱うので、都市群のデータの下のA28から貼り付けます。23人分のデータがA29からA51に入ります。
最初は、共通のパラメータを使うモデルを考えましょう。
農村部にもB列にlog(f(x_α|λ))の式を書いてください。共通のパラメータを使うのでλはC1セルの値を使います。
対数尤度はそれらの値を全部足したものですので、B52セルに =SUM(B2:B27)+SUM(B29:B51) と書いてください。

目的セル(E)欄に、最大にしたい対数尤度が書かれているB52、変化させるセル(B)欄にパラメータλが書かれているC1と書きます。欄をクリックしてからそれらのセルをクリックしても構いません。
今回は、パラメータλは正でなければなりませんので 制約のない変数を非負数にする(K) にチェックをつけたままにして解決(S)をクリックします。

うまく最適値が見つかりました。OKをクリックすると、対数尤度の最大値がB52に、そして対数尤度が最大になるλの値、つまり最尤推定値がC1に入力されます。

次に都市群と農村群で別のパラメータを用いるモデルを考えます。
E1セルに都市群のパラメータとしてとりあえず先ほど求めた共通のλの値を、E28セルに農村群のパラメータとしてとりあえず同じ値を入力して、D2セルにデータxとしては先ほど同様A2を、パラメータλとしてE1セルを使う式を書いてD27までコピーします。

今度はD29セルに、データxとしてはA29を、パラメータλとしてE28セルを使う式を書いてD51までコピーします。

D52セルにそれらの合計、つまり =SUM(D2:D27)+SUM(D29:D51) を書いて、この値を最大にするようにE1セルとE28セルの数字を変化させます。
SUM(D2:D27)を最大にするE1と、SUM(D29:D51)を最大にするE28を探せばよいのですが、これもソルバーを使いましょう。
最大にする目的セルはD52, 変化させるセルはE1,E28です。

対数尤度の最大値が求まりました。この大きさだけをみて、大きいほうが良い、と考えてはいけません。

AICを求めます。AIC=-2×対数尤度の最大値+2×パラメータ数です。これが小さいほうが良いモデルです。

別のパラメータを使うモデルのほうがAICが小さいので良いモデルといえます。実はこの場合は値の差が1以下なので断言し辛いですが、今回は練習ですので気にしないことにします。

ファイル操作

Excelで複数のファイルを操作する方法を説明しました。

一般化線形モデル

一般化線形モデルについて説明しました。

While文とIf文

For文と並んで重要な、While文とIf文について説明しました。

一般化情報量規準

一般化情報量規準を紹介しました。

Excel VBA入門

Excel VBAの導入を説明しました。
5月30日のExcel入門で使った身長と座高のデータを今回も使います。

生徒	身長	座高
1	159	88
2	150	84
3	157	86
4	153	81
5	158	83
6	152	85
7	155	83
8	157	83
9	145	76
10	158	85
11	161	85
12	150	83
13	148	79
14	154	84
15	154	85
16	159	85
17	149	83
18	155	86
19	153	84
20	160	88

出典

ヒストグラム

Excelを用いたヒストグラムの書き方を説明しました。

講義で使うデータ

120
106
112
168
120
140
160
160
116
118
104
172
130
124
134
142
112
112
110
120
134
136
138
144
120
120
160
108
140
98
142
130
138
126
108
120
136
118
118
122
128
124
110
134
116
138
116
158
106
166
118
158
124
124
110
166
132
126
122
132

出典

罰則付き最尤法

罰則付き最尤法の必要性について説明しました。

Excel入門

いつもの教室でExcelの使い方を説明しました。

来週も同じ教室です！
サンプルデータ

県名 男性 女性 滋賀県 699 716 京都府 1250 1359 大阪府 4256 4579 兵庫県 2645 2896 奈良県 649 726 和歌山県 457 515 出典:総務省統計局

生徒 身長 座高 1 159 88 2 150 84 3 157 86 4 153 81 5 158 83 6 152 85 7 155 83 8 157 83 9 145 76 10 158 85 11 161 85 12 150 83 13 148 79 14 154 84 15 154 85 16 159 85 17 149 83 18 155 86 19 153 84 20 160 88 出典

スプライン回帰

非線形回帰のひとつとしてスプライン回帰について説明しました。

気象データ

気象データの時系列解析について説明しました。

ソルバーを使って最尤推定量を求める

今日のテーマ

最尤推定量は(対数)尤度を最大にする値としてデータから計算できます。この計算式θ=θ(X₁,X₂,...,X_n)を具体的に書き表すことが出来る問題もあります(以前のレポートで出しました)が、計算式を求めることが出来ない問題もあります。
そこで今日はExcelを使って最尤推定量を求めて、AICを計算し、モデルを選ぶ方法を勉強しましょう。

準備

Excelを立ち上げたら左上の「ファイル」を押して、「オプション」をクリックします。そして左側の「アドイン」をクリックして、下の「管理:Excelアドイン」の横の「設定」をクリックします。

するとアドインの画面に切り替わるので「ソルバーアドイン」にチェックをつけてOKをクリックします。

サンプルデータ

男児 女児 期間(週) 体重(g) 期間(週) 体重(g) 40 2968 40 3317 38 2795 36 2729 40 3163 40 2935 35 2925 38 2754 36 2625 42 3210 37 2847 39 2817 41 3292 40 3126 40 3473 37 2539 37 2628 36 2412 38 3176 38 2991 40 3421 39 2875 38 2975 40 3231

出典:一般化線形モデル入門, Annette J.Dobson 著, 共立出版
これをコピーしてエクセルに貼り付けてください。散布図を書くとこのようになります。

考えるモデル

パラメータが少ない順番に

男女共通モデル

男女ともY=α+βX+ε, ε～N(0,σ²)
つまり妊娠期間Xに対して体重Yの確率分布は
男女ともN(α+βX,σ²)
推定するパラメータはα, β, σの3つ。

傾きだけ共通モデル

男児:Y=α₁+βX+ε, ε～N(0,σ²)
女児:Y=α₂+βX+ε, ε～N(0,σ²)
つまり妊娠期間Xに対して体重Yの確率分布は
男児:N(α₁+βX,σ²)
女児:N(α₂+βX,σ²)
推定するパラメータはα₁, α₂, β, σの4つ。

男女別モデル

男児:Y=α₁+β₁X+ε, ε～N(0,σ²)
女児:Y=α₂+β₂X+ε, ε～N(0,σ²)
つまり妊娠期間Xに対して体重Yの確率分布は
男児:N(α₁+β₁X,σ²)
女児:N(α₂+β₂X,σ²)
推定するパラメータはα₁, α₂, β₁, β₂, σの5つ。

まず最初に、パラメータが一番少ない男女共通モデルについて、α, β, σをExcelに求めてもらいます。
そのためには、まず最初に適当な初期値を与えます。するとExcelのソルバーはそこから出発して真の値に近づいていきます。初期値をどのように与えればよいでしょうか。最初から真の値が分かるならExcelのソルバーに頼る必要はないのですから、真の値に近くてしかも簡単に求められる値を初期値として使います。
直線y=α+βxのαとβの両方を簡単に知ることは出来ませんが、水平な線、つまりβ=0ならば、その高さαは平均値にしておけば大体データに近そうです。この時σは平均値からのずれ、つまり標準偏差にしておきましょう。
そこで早速体重の平均値と標準偏差を求めましょう。Excelに次のように貼り付けます。

A15に平均、A16に標準偏差と書いてから、B15に平均を求める式、B16に標準偏差を求める式を書きます。
B15にはB3からB14までとD3からD14までの平均ということで =AVERAGE(B3:B14,D3:D14) と書きます。先頭の等号を忘れないでください。
B16にはB3からB14までとD3からD14までの標準偏差ということで =STDEV(B3:B14,D3:D14) と書きます。
するとこのようになります。

セルの幅に合わせて四捨五入されますので、人によって表示される値が少し違うこともあります。
この値を参考にしてα, β, σの初期値を決めます。
まずA17, A18, A19セルにα, β, σと書いて
B17セルにはαの初期値として平均値を書きます。適当なところで四捨五入してかまいません。
B18セルにはβの初期値を書きます。傾きが良く分からなくて水平線を引くので0にします。
B19セルにはσの初期値として標準偏差を書きます。これも適当なところで四捨五入してかまいません。

対数尤度の計算

確率関数や密度関数にデータを代入した値が尤度、それにlogをつけたのが対数尤度です。独立な複数の観測値に関して尤度は掛け算になっているので、logをつけて掛け算を足し算にします。
log(Π_i=1ⁿf(X_i))=Σ_i=1ⁿlog(f(X_i))
正規分布N(μ,σ²)の密度関数は =NORMDIST(x, μ, σ, 0) です。最後の0を1に変えると分布関数になります。
今回の問題では、xのところが体重yで、平均値μがα+βxです。それにlogを付けますが、Excelでは自然対数は =LN(x) ですのでE3セルに =LN(NORMDIST(B3,B$17+B$18*A3,B$19,0)) と書きます。対応するセルに色がつくので確認してください。

α, β, σに相当するB17, B18, B19の数字の前には$を付けています。これはこの式をコピーして下に貼り付けても番号が変わらないようにするためです。
E4セルからE14セルにも同様に数式を書きます。一つ一つ書いていると手間がかかりますのでE3セルをコピーして貼り付けます。数式の中のセルの番号がどのように変わるか確認してください。
次にF3セルに女児の一人目の対数尤度の式 =LN(NORMDIST(D3,B$17+B$18*C3,B$19,0)) を書きます。

これもF4からF14へコピーして、E3からF14まで全部足したものが対数尤度です。
E15セルに対数尤度と書いて、F15セルに合計を求める式 =SUM(E3:F14) を書きます。

ソルバー

ここが今回のポイントです。これからα, β, σに色々な値を代入して対数尤度が最大になればよいのですが、試行錯誤するのは大変ですのでExcelのソルバーを使います。
エクセルの上の列の「データ」をクリックして、右側の「ソルバー」をクリックします。「ソルバー」が見当たらなかったら、このページの上の準備のところを見てください。

するとこのような画面が表示されるので、F15を最大にするためにB17からB19を変化させる、と指定します。
また、αは正の数とは限らないので「制約のない変数を非負数にする(K)」のチェックを外します。さらにソルバーの設定を変えるのでオプションを押します。

「GRG非線形」のタブを選んで、微分係数を「中央」にします。

OKを押してソルバーの画面に戻って「解決(S)」を押すと、対数尤度が大きくなるようにα, β, σの値が調整されて次のメッセージが表示されます。

OKをクリックすると、対数尤度の値が最大になっています。

AICの計算

このようにして求めた対数尤度の最大値からパラメーターの個数を引いた値が一番大きなモデルが、一番良いモデルなのですが、この方法が考えられた経緯から、-2倍して
-2×(対数尤度の最大値)+2×(パラメータの個数)
をAIC (Akaike information criterion)といい、マイナスを掛けているのでAICが一番小さなモデルが一番良いモデルです。
このモデルのAICは-2×(-159.265)+2×3=324.53です。

レポート課題

傾きだけ共通モデル、及び男女別モデルについて、最尤推定量及びAICの値を求め、どのモデルが一番良いか選びなさい。

赤池情報量規準

平均対数尤度を推定するために対数尤度を用いた時に生じるかたよりの計算を終わらせて、赤池情報量規準を紹介しました。

現時点での問題点

現時点で少々問題だと思える点を紹介しました。

検定

信頼区間と共に統計の基本的な考え方である、統計的仮説検定について説明しました。

偏りの評価

平均対数尤度を推定するために対数尤度を用いると、大きいほうに偏ってしまいます。どの程度偏るのかを計算します。

信頼区間

常識と理論の整合性をとるために、間違う確率をどの程度許せばどのような結論を導けるのかを説明しました。

kullback leibler情報量

最尤法とKullback-Leibler情報量の関係について説明しました。

ベイズの公式

最尤推定量

ポアソン分布と正規分布の最尤推定量の計算を説明しました。

データサイエンスへの招待

データサイエンスとはどのような学問なのか説明しました。

統計モデル

この講義の目的を、二つの例を用いて説明しました。