行列の演算

行列の演算について説明します。

線形代数への招待

高校の数学と大学の数学の違い、線形代数を学ぶ目的などを説明します。

多変量解析とは

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

出典