2013人文

2013 人文 (統計検定1級応用 過去問の解説です)

  1. 項目反応理論
  2. 多変量正規分布
  3. 主成分分析と因子分析
  4. 正方分割表の対称性の検定
  5. 重回帰分析
     

2013人文① 項目反応理論

項目反応理論を順序立って誘導解説した問題である。

教育心理学会 様 – 項目反応理論(IRT)の考え方と実践
     https://www.jstage.jst.go.jp/article/arepj/58/0/58_321/_pdf

(1) テストは n個の問題(項目)から構成される。
 各問題 \(i:(i=1,2,\cdots,n)\) に対し、確率変数 \(U_i\) は正解なら1、不正解なら 0 をとる。

各問題 \(i\) に対し、能力 \(\theta\) の生徒が正解する確率を、 \(P_i(\theta)=P(U_i=1|\theta)\) とすると、

不正解である確率は、\(Q_i(\theta)=1-P_i(\theta)\)

\(P_i(\theta)\) は項目特性関数といい、\(P\left( U_{i}=1| \theta \right) =c_i+\dfrac{1-c_i}{\left( 1+\exp \left[ -a_{i}\left(\theta -b_{i}\right) \right] \right) }\) で表される。

 \(b_i\) は困難度(項目に正答できるかできないか五分五分に必要な能力値)

 \(a_i\) は識別力(能力の違いをより敏感に識別する感度)

 \(b_i\) は、グラフの水平位置に対応し、 \(a_i\) はグラフの傾斜に対応する。

 \(c_i\) は、あて推量で正解する確率である。 

(2) 項目 \(i\) の配点を \(w_i\) とすると、テストの総得点Yは、

     \(Y=w_1U_1+\cdots +w_nU_n\)   で表される。

\(U_i\) はベルヌーイ分布に従い、

     \(E(U_i)=1\cdot P_i(\theta)+0\cdot Q_i(\theta)=P_i(\theta)\)

     \(V(U_i)=1^2\cdot P_i(\theta)-P_i(\theta)^2=P_i(\theta)Q_i(\theta)\)

よって、    \(\displaystyle E(Y)=\sum_{i=1}^n w_iP_i(\theta)\)    \(\displaystyle V(Y)=\sum_{i=1}^n w_i^2P_i(\theta)Q_i(\theta)\)

 

(3) 目的は能力 \(\theta\) を推定することである。

簡単のため、パラメータ \(a_i=a\:,\:\:\:b_i=b\:,\:\:\:c_i=c\) は既知で、全ての項目で共通の値とする。

仮に、被験者の正答数による正答率が \(z=\frac{k}{n}\) であるとすると、能力 \(\theta\) はどうなるだろうか?

パラメータが共通だから、項目特性関数も共通で、

     \(P_i(\theta)=P(\theta)\)    (すべての \(i\) で共通)

\(P(\theta)=\dfrac{k}{n}\) と予測され、 \(c+\dfrac{1-c}{\left( 1+\exp \left[ -a\left(\theta -b\right) \right] \right) }=\dfrac{k}{n}\) から \(\theta\) を求めることになる。

詳しくみていく。

局所独立性を仮定しているので、正答数は二項分布 \(B\left(n,P_i(\theta)\right)\) に従う。
(局所独立性とは、𝜃 を固定したとき各項目に対する解答が独立であることをいう)

能力 \(\theta\) のもとで正答数が \(k\) である確率は、

     \(L(\theta)={}_nC_kP(\theta)^k\left(1-P(\theta)\right)^{n-k}\)

これは尤度関数で、対数尤度関数は、

     \(l(\theta)=\log {}_nC_k+k\log P(\theta)+(n-k)\log \left(1-P(\theta)\right)\)

これが最大となるように \(\theta\) を決めるのが最尤法である。

     \(l'(\theta)=k\dfrac{P'(\theta)}{P(\theta)}+(n-k) \dfrac{-P'(\theta)}{1-P(\theta)} =0\) から、 \(P(\theta)=\dfrac{k}{n}=z\)   (予想どおり)

     \(c+\dfrac{1-c}{\left( 1+\exp \left[ -a\left(\theta -b\right) \right] \right) }=z\) から、\(\exp \left[ -a\left(\theta -b\right) \right]=\dfrac{1-z}{z-c}\)

     \(-a\left(\theta -b\right)=\log\dfrac{1-z}{z-c}\)     \(\therefore \left(\theta -b\right)=-\dfrac{1}{a}\log\dfrac{1-z}{z-c}=\dfrac{1}{a}\log\dfrac{z-c}{1-z}\)

推定値は    \(\widehat{\theta}=b+\dfrac{1}{a}\log\dfrac{z-c}{1-z}\)

(4) 一般的定義を確認する。

項目情報関数(Item Information Function)

     \(f_i(\theta)=a_i^2P_i(\theta)(1-P_i(\theta))\)   問題の情報量を表す。 \(P_i(\theta)=\dfrac{1}{2}\) で最大となる。

テスト情報関数(Test Information Function)

     \(\displaystyle I(\theta)=\sum_{i=1}^n f_i(\theta)\)     テスト全体の情報量を表す。

これが基本的定義であるが、本問では配点 としたテストの総得点

     \(Y=w_1U_1+\cdots +w_nU_n\)  のテスト情報関数を考える問題である。

誘導に従い、配点のある場合のテスト情報関数の定義を

     \(I\left( \theta \right) =\dfrac{\left( \dfrac{d}{d\theta }E\left(Y| \theta \right) \right) ^{2}}{V\left( Y| \theta \right) }\)     \(I(\theta)\) を最大にする配点は  \(w_{i}=\dfrac{P_{i}'\left( \theta \right) }{P_{i}\left( \theta \right) Q_{i}\left( \theta \right) }\) 

という事実を認めて話を進める。

あて推量で正解する確率 \(c_i=0\) のとき、 \(P_i(\theta) =\dfrac{1}{ 1+\exp \left[ -a_{i}\left(\theta -b_{i}\right) \right] }\) であるが、

興味深いことに、  \(P_i'(\theta)=a_iP_i(\theta)Q_i(\theta)\)  が成立する。

実際、 \(Q_i(\theta) =1-\dfrac{1}{ 1+\exp \left[ -a_{i}\left(\theta -b_{i}\right) \right] }=\dfrac{\exp \left[ -a_{i}\left(\theta -b_{i}\right) \right]}{ 1+\exp \left[ -a_{i}\left(\theta -b_{i}\right) \right] }\) より

     \(P_i'(\theta) =\dfrac{a_i\exp\left[ -a_{i}\left(\theta -b_{i}\right) \right] }{ \left(1+\exp \left[ -a_{i}\left(\theta -b_{i}\right) \right] \right)^2}=a_iP_i(\theta)Q_i(\theta)\)

よって、 \(\dfrac{P_{i}'\left( \theta \right) }{P_{i}\left( \theta \right) Q_{i}\left( \theta \right) }=a_i\)  これは \(\theta\) に依存しない定数である。

つまり、識別力を配点とすれば、テスト情報関数は最大となる。

 

2013人文② 多変量正規分布 

(1) 共分散は線形性を持つので、

\(\begin{split}
V\left(\dfrac{X_1+X_2+X_3}{3}\right)&=Cov\left(\dfrac{X_1+X_2+X_3}{3}\:,\:\dfrac{X_1+X_2+X_3}{3}\right)
=\dfrac{1}{9}Cov\left(X_1+X_2+X_3,X_1+X_2+X_3\right)\\
&=\dfrac{1}{9}\left(V(X_1)+V(X_2)+V(X_3)+2Cov(X_1,X_2)+2Cov(X_2,X_3)+2Cov(X_3,X_1)\right)\\
&=\dfrac{1}{9}\left(3\cdot 10^2+3\cdot 2\cdot 10^2\cdot0.7\right)=80
\end{split}\)

あとは解答を参照。

(2)

[多変量正規分布の基本事項]

n次元標準正規分布 \(Z\sim N(0,I_n)\) のとき、正則行列 Aとμを用いて、\(X=AZ+\mu\) とすると、

     \(X\sim N(\mu,AA')\)     \(AA'=\Sigma\) は対称な正定値行列

     確率密度関数         \(f(x)=\dfrac{1}{(2\pi)^n\sqrt{|\Sigma|}}\exp\left({-\dfrac{1}{2}X'\ \Sigma^{-1}X}\right)\)

逆に、Xの各成分が1次元正規分布に従い、分散共分散行列Σが正則であるとき、

     Σは対称な正定値行列であり、 \(AA'=\Sigma\) となる正則行列Aが存在。

     \(Z=A^{-1}(X-\mu)\) により、n次元標準正規分布Zに変換される。

[コレスキー分解]

一般に、正定値対称行列Aは、適当な下三角行列L(上三角行列でも可)を用いて

     \(A=LL'\)     のように分解できる。(コレスキー分解という)

3次元の場合で成分計算する。

\({A =
\begin{pmatrix}
a_{11} & a_{21} & a_{31} \\
a_{21} & a_{22} & a_{32} \\
a_{31} & a_{32} & a_{33} \\
\end{pmatrix},\quad
L =
\begin{pmatrix}
\, l_{11} & 0 & 0 \\
\, l_{21} & l_{22} &0\\
\, l_{31} & l_{32} & l_{33} \\
\end{pmatrix}
}\) とおくと、\(
LL^{\rm T} =
\begin{pmatrix}
l_{11} & 0 & 0 \\
l_{21} & l_{22} & 0 \\
l_{31} & l_{32} & l_{33} \\
\end{pmatrix}

\begin{pmatrix}
l_{11} & l_{21} &\: l_{31} \\
0 & l_{22} &\:l_{32} \\
0 & 0 &l_{33} \\
\end{pmatrix}
\)

1列目を比べて、 \(a_{11} = l_{11}^2\;,\:a_{21} = l_{11}l_{21}\:,\: a_{31} = l_{11}l_{31}\)

これより \(l_{11}=\sqrt{a_{11}}\:\:,\:\:l_{21}=\dfrac{a_{21}}{l_{11}}\:\:,\:\:l_{31}=\dfrac{a_{31}}{l_{11}}\) が順に求まる。

\(a_{22} = l_{21}^2+l_{22}^2\) より 、 \(l_{22}=\sqrt{a_{22}-l_{21}^2}\) 、   \(a_{32} = l_{21}l_{31}+l_{22}l_{32}\) より \(l_{32}=\dfrac{a_{32}- l_{21}l_{31}}{l_{22}}\)

\(a_{33} = l_{31}^2+l_{32}^2+l_{33}^2
\)  より \(l_{33}=\sqrt{a_{33}-l_{31}^2-l_{32}^2}\)   順に求まる。

また、正定値対称行列Aについて、次の4つは同値である。

   (a) 0 でない任意の \(x\in R^n\) に対して \(x'Ax>0\) が成り立つ。

   (b) \(A\) の固有値はすべて正である。

   (c) ある直交行列 \(Q\) と対角要素がすべて正の対角行列 \(D\) を用いて \(A=QDQ'\) と書ける。

   (d) \(A=SS'\) を満たす正則行列 \(S\) が存在する。

詳しくは、東京大学 寒野 善博 様
https://www.or.mist.i.u-tokyo.ac.jp/kanno/lecture/math_design/psd170420.pdf  参照

本問では、

\((X_1,X_2,X_3)\) は3変量正規分布に従っており、分散共分散行列 \(\Sigma\) は正定値対称行列である。

よって、\(\Sigma=AA'\) となる正則行列 \(A\) が存在し、コレスキー分解により求めることができる。

\(l_{11}=\sqrt{a_{11}}=10\:\:,\:\:l_{21}=\dfrac{a_{21}}{l_{11}}=\dfrac{10^2\cdot0.7}{10}=7\:\:,\:\:l_{31}=\dfrac{a_{31}}{l_{11}}=7\)

\(l_{22}=\sqrt{a_{22}-l_{21}^2}=\sqrt{10^2-7^2}=\sqrt{51}\)   \(l_{32}=\dfrac{a_{32}- l_{21}l_{31}}{l_{22}}=\dfrac{10^2\cdot0.7- 7^2}{\sqrt{51}}=\dfrac{21}{\sqrt{51}}\)

\(l_{33}=\sqrt{a_{33}-l_{31}^2-l_{32}^2}=\sqrt{10^2-7^2-\dfrac{21^2}{51}}=\sqrt{51-\dfrac{21^2}{51}}=\sqrt{\dfrac{51^2-21^2}{51}}=\sqrt{\dfrac{2650}{51}}\)

\(A=\begin{pmatrix}
10&0&0\\
7&\sqrt{51}&0\\
7&\dfrac{21}{\sqrt{51}}&\dfrac{\sqrt{2160}}{\sqrt{51}}
\end{pmatrix}\) となり、  \(Z=A^{-1}(X-\mu)\) より \(X=AZ+\mu\) だから

独立な標準正規乱数 \((Z_1,Z_2,Z_3)\) に対し、

     \(X_1=10Z_1+50\) 、\(X_2=7Z_2+\sqrt{51}Z_2+50\) 、\(X_3=7Z_1+\dfrac{21}{\sqrt{51}}Z_2+\dfrac{\sqrt{2160}}{\sqrt{51}}Z_3\)

とすればよい。(これは公式問題集の解答と一致する)

公式問題集では4つの標準正規乱数を用いる方法が紹介されているが、コレスキー分解のほうが一般の場合に通用する方法である。

[参考サイト]   

関口 叡範 様 「Cholesky分解を用いた正定値対称行列の効率的な逆行列の計算方法」
     https://qiita.com/sekky0816/items/6f019263891c6e0b36d0

千葉大学工学部 武居 昌宏 様 「数値計算法」
     http://www.em.eng.chiba-u.jp/~takei/NOTE/suutikeisan/4-5-6_20170418.pdf

(3) 合格基準 C の合格率が 0.4 であるとする。

n回試行したときの合格回数を M とすると、\(M\sim B(n,0.4)\) に従い、

正規分布で近似すると、\(M\sim N\left(0.4\cdot n\:\:,\:\:(0.4)(0.6)n\right)\)

合格率は近似的に \(\dfrac{M}{n}\sim N\left(0.4\:\:,\:\:\dfrac{(0.4)(0.6)}{n}\right)\)

観測値 \(\dfrac{M}{n}\) は95%の確率で、 \(\left(0.4-1.96\sqrt{\dfrac{0.24}{n}}\:\:,\:\:0.4+1.96\sqrt{\dfrac{0.24}{n}}\right)\) に入るはずである。

判断の精度を上げるには信頼区間の幅を小さくする必要がある。

標準誤差は予測値の標準偏差であり、

     \(\sqrt{\dfrac{0.24}{n}} \leq0.01 \) とすると \(n\geq2400\)

よって、2400回以上行えば、標準誤差は0.01に抑えられる。

 

2013人文③ 主成分分析と因子分析

表2について説明する。

変数を標準化した場合
 分散共分散行列は相関行列になる。
 固有値は新変数の分散であり、分散合計は変換後も不変であるから、
 (固有値の合計)=(変数の個数)で 6 になっている。

変数を標準化していない場合
 分散共分散行列を対角化する直交行列から新変量がつくられる。
  各変数の分散は大きいままだから、固有値も大きいままである。
 相関行列を用いても新変量は構成できるが、新変量の分散共分散行列が対角行列にならない。
  標準化されていない場合は分散共分散行列を使わなければならない。

どちらの場合でも、次のことが成立している。

・固有ベクトルが主成分への変換式の係数であり、固有値は新変数の分散である。
・新しい変量の分散共分散行列は固有値を対角成分にもつ対角行列である。
・新しい変数は無相関であり、第1主成分は新しい変量で分散が最大のものである。

(1) 変数が標準化されていない分散共分散行列を用いた主成分分析である。

直交行列による対角化なので、固有ベクトルの表の列の平方和は1である。

第1主成分は、 \(F_1=(0.398)PC1+\cdots+(0.441)RW3\)

分散は固有値に等しいから、\(V(F_1)=948.19\)

変数の平均は基準化されているので、\(E(F_1)=0\)

正規分布の和は正規分布に従うので、\(F_1\sim N(0\:,\:948.19)\)

     \(Z=\dfrac{F_1}{\sqrt{948.19}}\sim N(0,1)\) であるので、上側20%点を正規分布表から求めると、

     \(F_1=0.84\times \sqrt{948.19}\) と計算される。

(2) 公式問題集ではデータ行列を用いて説明しているが、ここでは確率変数を用いて説明する。

PC1からRW3までの6つのデータを成分とする6行1列の行列を \(X\) とする。

基準化され平均は0であるので、分散共分散行列 \(S\) は

     \(S=E(XX')\) である。( \(E\) は期待値を意味する)

\(S\) の最大固有値 \(\lambda_1\) に対応する固有ベクトルを \(h\) 、すべての成分が1のベクトルを \(l\) とするとき、

     第1主成分は \(F_1=h'X\)      6科目の合計点は  \(w=l'X\)

基準化されているので、

     \(Cov(w,F_1)=E(w\cdot F_1')=E(l'XX'h)=l'E(XX')h=l'Sh=l'\lambda_1 h=\lambda_1 l'h\)

\(l'h\) は固有ベクトルの成分の和であるから、\(l'h=2.406\)

固有値 \(\lambda_1\) は \(F_1\) の分散に等しいから \(s(F_1)=\sqrt{\lambda_1}\)

表1より、合計点 \(w\) の標準偏差は \(s(w)=74.31\)

よって、\(r(w,F_1)=\dfrac{Cov(w,F_1)}{s(w)s(F_1)}=\dfrac{\lambda_1 l'h}{s(w)\sqrt{\lambda_1}}=\dfrac{\sqrt{\lambda_1}l'h}{s(w)}=\dfrac{\sqrt{948.19}\cdot (2.406)}{74.31}=0.9970\)

相関が高く、直線関係にある。よって、合計点で代用することができる。

 

(3) 主成分分析を用いて因子分析に進めてみる。

相関行列を用いた、固有値と固有ベクトルの表を見る。

相関行列は直交行列で対角化できる。

     固有値と固有ベクトル

1 2 3 4 5 6 合計
固有値 3.926 1.038 0.352 0.321 0.190 0.172 6.0
固有ベクトル 平方和
0.438 -0.334 0.085 0.212 0.399 -0.696 1.0
0.416 -0.425 -0.227 -0.111 -0.762 -0.033 1.0
0.430 -0.386 0.093 0.091 0.378 0.711 1.0
0.384 0.448 0.687 0.307 -0.291 0.037 1.0
0.366 0.503 -0.679 0.380 0.064 0.058 1.0
0.411 0.323 0.005 -0.834 0.166 -0.055 1.0
平方和 1.0 1.0 1.0 1.0 1.0 1.0

直交行列なので各行、各列の平方和は1になっている。

固有値はそれぞれの主成分の分散であり、合計は元の変数の個数に一致する。

因子分析(主因子法)では、各因子の分散を1にしなければならない。

各列に固有値の平方根をかけたものが因子負荷行列になる。

(当サイト「負荷量と寄与率」を参照)

     因子負荷行列

固有値の平方根
1.981 1.019 0.593 0.567 0.436 0.415
平方和
因子負荷行列 0.868 -0.340 0.050 0.120 0.174 -0.289 1.0
0.824 -0.433 -0.135 -0.063 -0.332 -0.014 1.0
0.852 -0.393 0.055 0.052 0.165 0.295 1.0
0.761 0.456 0.408 0.174 -0.127 0.015 1.0
0.725 0.512 -0.403 0.215 0.028 0.024 1.0
0.814 0.329 0.003 -0.473 0.072 -0.023 1.0
平方和 3.927 1.037 0.352 0.321 0.190 0.172
累積 3.927 4.964 5.316 5.637 5.827 5.998
累積寄与率 0.654 0.827 0.886 0.939 0.971 1.000

各列の平方和は固有値になっている。

元の変数と因子は標準化されているので、各行の平方和は1となる。

第2因子で累積寄与率が80%を超えるので、因子数は2とする。

詳細は解答参照。

(4) 上の表は主因子法による。主因子法では、元の変数は無相関な変数(直交する変数)に変換される。固有値は影響力の大きさを表すので、固有値を降順に並べ、寄与率の累積が8割を越えれば、その個数の変数で十分説明できていると考える。

因子負荷行列の成分は、その変数と因子との相関係数であるから、相関係数の小さいものは排除することが望ましい。そこで因子負荷行列が単純構造に近くなるように回転を行う。

回転については、当サイト「因子分析は未定係数法」を参照。

 

2013人文④ 正方分割表の対称性の検定 

分割表の独立性の検定はよく知られている。(当サイト「分割表の分析」参照)

独立でないことが分かった段階で、対称性があるか検定したい場合がある。

本問では対象 \((i,j)\) の相互感情は独立でないことは予想できるが、対称性があるかは疑問である。

つまり、自分が相手に好意を持てば相手も好意を持ってくれるかという問題である。

解説するサイトはあまりないようであるが、次のサイトが参考になるかもしれない。

東京理科大学 田畑耕治様・富澤貞男様  「多元分割表統計解析における対称性のモデリング」
     https://www.jstage.jst.go.jp/article/sugaku/63/3/63_0633338/_pdf/-char/ja

公式問題集では丁寧な誘導と解説があるが、ここでは簡単のため次元を下げて考えてみたい。

イメージを持つために、日米の相互感情の分割表を考えることにしよう。

\(M=2\) 即ち、感情は(好き、嫌い)の2値とする。

\((f_{121},f_{122})\)
\((f_{211},f_{212})\)

(1) 分割表の \((i,j)\) として、

\((1,2)\) は日本が米国に対して持つ感情の調査であり、

     調査人数: \(n_{12}\)    (好き,嫌い)の度数: \((f_{121},f_{122})\)    (好き,嫌い)の確率: \((p_{121},p_{122})\)

            ( \(f_{121}+f_{122}=n_{12}\:\:,\:\:p_{121}+p_{122}=1\) )

     \((f_{121},f_{122})\) の尤度は  \(\displaystyle \dfrac{n_{12}}{f_{121}!\cdot f_{122}!}p_{121}^{f_{121}}\cdot p_{122}^{f_{122}}\)   これは二項分布である。

\((2,1)\) は米国が日本に対して持つ感情の調査であり、

     調査人数: \(n_{21}\)    (好き,嫌い)の度数:\((f_{211},f_{212})\)    (好き,嫌い)の確率:\((p_{211},p_{212})\)

            ( \(f_{211}+f_{212}=n_{21}\:\:,\:\:p_{211}+p_{212}=1\) )

     \((f_{211},f_{212})\) の尤度は  \(\displaystyle \dfrac{n_{21}}{f_{211}!\cdot f_{212}!}p_{211}^{f_{211}}\cdot p_{212}^{f_{212}}\)   これも二項分布である。

同時尤度は、

     \(\displaystyle L=\left(\dfrac{n_{12}}{f_{121}!\cdot f_{122}!}p_{121}^{f_{121}}\cdot p_{122}^{f_{122}}\right)\cdot \left(\dfrac{n_{21}}{f_{211}!\cdot f_{212}!}p_{211}^{f_{211}}\cdot p_{212}^{f_{212}}\right)\)

対数尤度は

\(\displaystyle \begin{split}\log L&=\log\left(\dfrac{n_{12}}{f_{121}!\cdot f_{122}!}\dfrac{n_{21}}{f_{211}!\cdot f_{212}!}\right)+\log \left(p_{121}^{f_{121}}\cdot p_{122}^{f_{122}}\cdot p_{211}^{f_{211}}\cdot p_{212}^{f_{212}}\right)\\
&=\log\left(\dfrac{n_{12}}{f_{121}!\cdot f_{122}!}\dfrac{n_{21}}{f_{211}!\cdot f_{212}!}\right)
+{f_{121}}\log \left(p_{121}\right)+{f_{122}}\log \left(p_{122}\right)+{f_{211}}\log \left(p_{211}\right)+{f_{212}}\log \left(p_{212}\right)

\end{split}\)

\(p_{121}+p_{122}=1\:\:,\:\:p_{211}+p_{212}=1\) のもとでの対数尤度の最大を考えるが、各 \((i,j)\) の調査は独立なので、

     \(p_{121}+p_{122}=1\) のもとでの \({f_{121}}\log \left(p_{121}\right)+{f_{122}}\log \left(p_{122}\right)
\)

     \(p_{211}+p_{212}=1\) のもとでの \({f_{211}}\log \left(p_{211}\right)+{f_{212}}\log \left(p_{212}\right)
\)

別々に最大を考えればよいことになる。

ラグランジュの未定乗数法により、

     \(Q_{12}={f_{121}}\log \left(p_{121}\right)+{f_{122}}\log \left(p_{122}\right)
-\lambda_{12}(p_{121}+p_{122}-1)\) を \(p_{121}\) で偏微分すると

     \(\dfrac{f_{121}}{p_{121}}-\lambda_{12}=0\)   \(\therefore f_{121}=\lambda_{12}\cdot p_{121}\)    同様に、\(f_{122}=\lambda_{12}\cdot p_{122}\)

\(n_{12}=f_{121}+f_{122}=\lambda_{12}\cdot p_{121}+\lambda_{12}\cdot p_{122}=\lambda_{12}(p_{121}+p_{122})=\lambda_{12}\) となるので、

     \(p_{121}=\dfrac{f_{121}}{n_{12}}\:\:,\:\:p_{122}=\dfrac{f_{122}}{n_{12}}\) が最尤推定量となる。

同様にして、 \(p_{211}=\dfrac{f_{211}}{n_{21}}\:\:,\:\:p_{212}=\dfrac{f_{212}}{n_{21}}\) が得られる。これは自然な結果である。

最大対数尤度は、

\(\log L_{max}=
\log\left(\dfrac{n_{12}}{f_{121}!\cdot f_{122}!}\dfrac{n_{21}}{f_{211}!\cdot f_{212}!}\right)
+{f_{121}}\log \left(\dfrac{f_{121}}{n_{12}}\right)+{f_{122}}\log \left(\dfrac{f_{122}}{n_{12}}\right)+{f_{211}}\log \left(\dfrac{f_{211}}{n_{21}}\right)+{f_{212}}\log \left(\dfrac{f_{212}}{n_{21}}\right)
\)

(2) 帰無仮説:\((p_{121},p_{122})=(p_{211},p_{212})\) のもとでは、

\(\log L^{(0)}
=\log\left(\dfrac{n_{12}}{f_{121}!\cdot f_{122}!}\dfrac{n_{21}}{f_{211}!\cdot f_{212}!}\right)
+({f_{121}}+f_{211})\log \left(p_{121}\right)+({f_{122}}+{f_{212}})\log \left(p_{122}\right)\)

ラグランジュの未定乗数法により最尤推定を行う。

\(Q_{12}=({f_{121}}+f_{211})\log \left(p_{121}\right)+({f_{122}}+{f_{212}})\log \left(p_{122}\right)-\lambda_{12}(p_{121}+p_{122}-1)\)  を \(p_{121}\) で偏微分すると

     \(\dfrac{{f_{121}}+f_{211}}{p_{121}}-\lambda_{12}=0\)  \(\therefore {f_{121}}+f_{211}=\lambda_{12}\cdot p_{121}\)  同様に \(f_{122}+f_{212}=\lambda_{12}\cdot p_{122}\)

\(\begin{split}
n_{12}+n_{21}&=(f_{121}+f_{122})+(f_{211}+f_{212})\\
&=(f_{121}+f_{211})+(f_{122}+f_{212})=\lambda_{12}\cdot p_{121}+\lambda_{12}\cdot p_{122}=\lambda_{12}(p_{121}+p_{122})=\lambda_{12}
\end{split}\)

よって、\(p_{121}=\dfrac{f_{121}+f_{211}}{n_{12}+n_{21}}\:\:,\:\:p_{122}=\dfrac{f_{122}+f_{212}}{n_{12}+n_{21}}\)

帰無仮説のもとでの最大対数尤度は、

\(\log L_{max}^{(0)}=
\log\left(\dfrac{n_{12}}{f_{121}!\cdot f_{122}!}\dfrac{n_{21}}{f_{211}!\cdot f_{212}!}\right)
+({f_{121}}+f_{211})\log \left(\dfrac{f_{121}+f_{211}}{n_{12}+n_{21}}\right)+({f_{122}}+{f_{212}})\log \left(\dfrac{f_{122}+f_{212}}{n_{12}+n_{21}}\right)\)

(4) 尤度比検定の考え方については

関西大学 浅野晃 先生 「検出力と尤度比検定」
          http://racco.mikeneko.jp/Kougi/2011a/AAN/2011aaan14.pdf 参照

尤度比は \(\Lambda=\dfrac{ L_{max}^{(0)}}{L_{max}}\) であるから、

     \(\begin{split}\log\Lambda&=\log{ L_{max}^{(0)}}-{\log L_{max}}\\
&=({f_{121}}+f_{211})\log \left(\dfrac{f_{121}+f_{211}}{n_{12}+n_{21}}\right)+({f_{122}}+{f_{212}})\log \left(\dfrac{f_{122}+f_{212}}{n_{12}+n_{21}}\right)\\
&\:\:\:\:-{f_{121}}\log \left(\dfrac{f_{121}}{n_{12}}\right)-{f_{122}}\log \left(\dfrac{f_{122}}{n_{12}}\right)-{f_{211}}\log \left(\dfrac{f_{211}}{n_{21}}\right)-{f_{212}}\log \left(\dfrac{f_{212}}{n_{21}}\right)\\
&={f_{121}}\left[\log \left(\dfrac{f_{121}+f_{211}}{n_{12}+n_{21}}\right)-\log \left(\dfrac{f_{121}}{n_{12}}\right)\right]
+{f_{122}}\left[\log \left(\dfrac{f_{122}+f_{212}}{n_{12}+n_{21}}\right)-\log \left(\dfrac{f_{122}}{n_{12}}\right)\right]\\
&\:\:\:\:+{f_{211}}\left[\log \left(\dfrac{f_{121}+f_{211}}{n_{12}+n_{21}}\right)-\log \left(\dfrac{f_{211}}{n_{21}}\right)\right]
+{f_{212}}\left[\log \left(\dfrac{f_{122}+f_{212}}{n_{12}+n_{21}}\right)-\log \left(\dfrac{f_{212}}{n_{21}}\right)\right]\\
&={f_{121}}\left[\log \left(\widehat{p}_{121}^{(0)}\right)-\log \left(\dfrac{f_{121}}{n_{12}}\right)\right]
+{f_{122}}\left[\log \left(\widehat{p}_{122}^{(0)}\right)-\log \left(\dfrac{f_{122}}{n_{12}}\right)\right]\\
&\:\:\:\:+{f_{211}}\left[\log \left(\widehat{p}_{211}^{(0)}\right)-\log \left(\dfrac{f_{211}}{n_{21}}\right)\right]
+{f_{212}}\left[\log \left(\widehat{p}_{212}^{(0)}\right)-\log \left(\dfrac{f_{212}}{n_{21}}\right)\right]\\
&={f_{121}}\left[\log \left(n_{12}\widehat{p}_{121}^{(0)}\right)-\log \left({f_{121}}\right)\right]
+{f_{122}}\left[\log \left(n_{12}\widehat{p}_{122}^{(0)}\right)-\log \left({f_{122}}\right)\right]\\
&\:\:\:\:+{f_{211}}\left[\log \left(n_{21}\widehat{p}_{211}^{(0)}\right)-\log \left({f_{211}}\right)\right]
+{f_{212}}\left[\log \left(n_{21}\widehat{p}_{212}^{(0)}\right)-\log \left({f_{212}}\right)\right]
\end{split}
\)

よって、

     \(\begin{split}-\log\Lambda=&{f_{121}}\left[\log \left({f_{121}}\right)-\log \left(n_{12}\widehat{p}_{121}^{(0)}\right)\right]
+{f_{122}}\left[\log \left({f_{122}}\right)-\log \left(n_{12}\widehat{p}_{122}^{(0)}\right)\right]\\
&+{f_{211}}\left[\log \left({f_{211}}\right)-\log \left(n_{21}\widehat{p}_{211}^{(0)}\right)\right]
+{f_{212}}\left[\log \left({f_{212}}\right)-\log \left(n_{21}\widehat{p}_{212}^{(0)}\right)\right]
\end{split}\)

分割表での対数尤度比統計量は \(-2\log\Lambda\) であり、漸近的にカイ二乗分布に従う。

Wikipedia様 https://ja.wikipedia.org/wiki/尤度比検定 参照)

自由度は、事前に決めたパラメーターの数を減ずることに注意する必要がある。

統計科学研究所様 「適合度検定」(p5)参照)
         https://statistics.co.jp/reference/StatisticalDataScienceA/Chapter7nonpassword.pdf
未知母数を含む分布の検定する場合、推定した未知母数の数だけ自由度が減る。

尤度比検定については統計数理⑤でも出題されている。

今後、さらに調べて補足する予定であるが、初見で解答するのは難しい問題である。

 

2013人文⑤ 重回帰分析

(1)は単回帰分析による分析

単回帰モデル \(y=\alpha+ \beta x+u\) のOLS推定量は、

     \(\displaystyle \beta=\frac{\sum ^{n}_{i=1}\left( x_{i}-\overline{x}\right) \left( y_{i}-\overline{y}\right) }{
\sum ^{n}_{i=1}\left( x_{i}-\overline{x}\right) ^{2}}=\dfrac{s_{xy}}{s_x^2}=\dfrac{s_{xy}}{s_xs_y}\cdot\dfrac{s_y}{s_x}=r_{xy}\cdot\dfrac{s_y}{s_x}\)

決定係数は、

     \(\begin{split}
R^2&=1-\frac
{\sum ^{n}_{i=1}\left( y_{i}-(\alpha+ \beta x_i)\right)^2 }
{\sum ^{n}_{i=1}\left( y_{i}-\overline{y}\right) ^{2}}
=\frac
{\sum ^{n}_{i=1}\left( (\alpha+ \beta x_i)-\overline{y}\right)^2 }
{\sum ^{n}_{i=1}\left( y_{i}-\overline{y}\right) ^{2}}\\
&=\frac
{\sum ^{n}_{i=1}\left( (\alpha+ \beta x_i)-\overline{y}\right)^2 }
{\sum ^{n}_{i=1}\left( y_{i}-\overline{y}\right) ^{2}}
=\frac
{\sum ^{n}_{i=1}\left( (\alpha+ \beta x_i)-(\alpha+\beta\overline{x})\right)^2 }
{\sum ^{n}_{i=1}\left( y_{i}-\overline{y}\right) ^{2}}\\
&=\beta^2\frac
{\sum ^{n}_{i=1}\left( x_{i}-\overline{x}\right) ^{2}}
{\sum ^{n}_{i=1}\left( y_{i}-\overline{y}\right) ^{2}}
=r_{xy}^2\cdot\dfrac{s_y^2}{s_x^2}\cdot\dfrac{s_x^2}{s_y^2}
=r_{xy}^2
\end{split}\)

     [参照]高校数学の美しい物語様 https://manabitimes.jp/math/1016

単回帰式 \(y=a_0+a_1x_1\) を推定すると、

     \(a_1=1.294\)  、\(a_0=-2.359\)  、\(y=-2.359+1.294x_1\)  、\(R^2=(0.772)^2=0.596\)

単回帰式 \(y=c_0+c_1x_2\) を推定すると、

     \(c_1=0.467\)  、\(c_0=34.837\)  、\(y=34.837+0.467x_2\)  、 \(R^2=(0.344)^2=0.118\) 

(2)は単回帰と重回帰をモデルとして比較する。

単回帰の自由度調整済み決定係数

     \(\begin{split}R^{*2}&=1-\frac
{\sum ^{n}_{i=1}\left( y_{i}-(\alpha+ \beta x_i)\right)^2 /(n-k-1)}
{\sum ^{n}_{i=1}\left( y_{i}-\overline{y}\right) ^{2}/(n-1)}\\
&=1-\dfrac{S_e/(n-k-1)}{S_{yy}/(n-1)}
=1-\left(\dfrac{S_e}{S_{yy}}\right)\left(\dfrac{n-1}{n-k-1}\right)
\end{split}\)

     ( \(k\) は説明変数の個数、単回帰では \(k=1\) )

全変動 \(\displaystyle S_{yy}=\sum ^{n}_{i=1}\left( y_{i}-\overline{y}\right) ^{2}\)      回帰変動 \(\displaystyle S_R=\sum ^{n}_{i=1}\left( (\alpha+ \beta x_i)-\overline{y}\right)^2\)

残差変動  \(\displaystyle S_e=\sum ^{n}_{i=1}\left( y_{i}-(\alpha+ \beta x_i)\right)^2 \) であり、\(S_{yy}=S_R+S_e\) 、\(R^2=\dfrac{S_R}{S_{yy}}\) が成り立つ。

単回帰については分析表が与えられていないので、相関係数などを利用して計算する。

\(S_e=S_{yy}-S_R=S_{yy}-R^2\cdot S_{yy}=(1-R^2)S_{yy}\)  から、 \(\dfrac{S_e}{S_{yy}}=1-R^2\)  を利用する。

単回帰(2)について、 \(R^2=(0.772)^2=0.596\) より

     \(R^{*2}=1-\left(1-R^2\right)\left(\dfrac{n-1}{n-k-1}\right)
=1-(1-0.596)\cdot \dfrac{119}{118}=0.593\)

\(0.593<0.620\)  であるので、重回帰分析の方が当てはまりが良い。

(3)解答参照

[更新日]

© hyoccori-toukei 2021