人文2015

2015 人文　（統計検定1級応用　過去問の解説です）

項目反応理論
主成分分析と因子分析
多次元尺度構成法（MDS）
ベイズ予測
２元配置分散分析

2015人文① 項目反応理論

パラメーターとしては、

識別力（a：能力の検出力）、困難度（b：問題の難度）がある。

a はグラフの傾きを、b はグラフの水平位置を表す。

能力 \(\theta\) の生徒が問題 \(j\) に正解する確率を \(p_j(\theta)\) とすると、

　・プロビットモデル（正規分布の累積関数を使用）　　　\(p_j(\theta)=\Phi\left(a_j(\theta-b_j)\right)\)　

　・ロジスティックモデル（指数関数を使用）　　　　　　　　\(p_j(\theta) =\dfrac{1}{ 1+\exp \left[ -Da_{j}\left(\theta -b_{j}\right) \right]}\)

確率変数 \(x_j=0,1\)　：項目 \(j\) に対する反応結果を表し、 \(x_j\) は独立とする。

１母数モデル

　　　　　識別力 \(a_j=a\) （一定）：グラフの傾きが一定のモデル

２母数モデル

　　　　　識別力 \(a_j\) 、困難度　\(b_j\) が j に依存するモデル

(1)　１母数プロビットモデルで、\(a=1,\ b=0\) のとき

m問中k問正解したときの能力θについて最尤推定量を求める。

尤度関数は　　　　　\(L(\theta)={}_mC_k \Phi(\theta)^k\left(1-\Phi(\theta)\right)^{m-k}\)

対数尤度関数は　　\(l(\theta)=\log {}_mC_k+k\log \Phi(\theta)+(m-k)\log \left(1-\Phi(\theta)\right)\)

\(\Phi'(\theta)=\phi(\theta)\) として、\(l'(\theta)=k\dfrac{\phi(\theta)}{\Phi(\theta)} -(m-k)\dfrac{\phi(\theta)}{1-\Phi(\theta)} =0\)

\(k(1-\Phi(\theta))=(m-k)\Phi(\theta)\)　　　　　\(\therefore k=m\Phi(\theta)\)

よって、 \(\Phi(\theta)=\dfrac{k}{m}\) （ \(\Phi(\theta)\) は正解する確率だから、妥当な結果といえる）

　　　　　\(\widehat\theta=\Phi^{-1}\left(\dfrac{k}{m}\right)\) 分布関数のグラフで \(\dfrac{k}{m}\) に対応する値が最尤推定値である。

(2) １母数ロジスティックモデルで、 m問中k問正解したときの能力θの最尤推定量を求める。

\(a=1,\ b=0\) 、\(p(\theta) =\dfrac{1}{ 1+\exp \left[ -D\theta \right]}\) 、尤度関数は \(L(\theta)={}_mC_k p(\theta)^k\left(1-p(\theta)\right)^{m-k}\)

対数尤度関数は

\(\begin{split}
l(\theta)&=\log {}_mC_k+k\log\left[\dfrac{1}{ 1+\exp \left[ -D\theta \right]}\right]+(m-k)\log\left[ \dfrac{\exp \left[ -D\theta \right]}{ 1+\exp \left[ -D\theta \right]}\right]\\
&=\log {}_mC_k-k\log\left[{ 1+\exp \left[ -D\theta \right]}\right]+(m-k) \log \left[\exp \left[ -D\theta \right] \right]
+(-m+k)\log\left[ { 1+\exp \left[ -D\theta \right]}\right]\\
&=\log {}_mC_k-m \log \left[1+\exp \left[ -D\theta \right] \right]-(m-k)D\theta
\end{split}\)

\(\therefore l'(\theta)=m\dfrac{D\exp \left[ -D\theta \right]}{1+\exp \left[ -D\theta \right]} -(m-k)D=0\) 　　　　　より、\(m\left(1-\dfrac{1}{1+\exp \left[ -D\theta \right]}\right )=
m-k\)

よって、(\dfrac{1}{1+\exp \left[ -D\theta \right]}=\frac{k}{m}\)　　　　　（これは \(p(\theta)=\dfrac{k}{m}\) を意味するから妥当といえる）

\(\theta\) の最尤推定値を \(\widehat \theta\) とすれば、 \(1+\exp(-D\widehat\theta)=\dfrac{m}{k}\)　　　　　\(\therefore \exp(-D\widehat\theta)=\dfrac{m-k}{k}\)

　　　　　\(-D\widehat\theta=\log \left(\dfrac{m-k}{k}\right)\) 　　　　よって、　 \(\widehat\theta=\dfrac{1}{D}\log \left(\dfrac{k}{m-k}\right)\)

(3) 実際のデータが与えられたとき、能力 \(\theta\) の推定値を求める。

正答数が8のとき、\(m=10~,~k=8\)

プロビットモデルでは、

　　　　　\(\widehat\theta=\Phi^{-1}\left(\dfrac{k}{m}\right)=\Phi^{-1}\left(0.8\right)\)

　　　　　正規分布表で、上側確率 0.2 となる点を調べると、\(\hat{\theta}=0.84\)

ロジスティックモデルでは、

　　　　　\(\widehat\theta=\dfrac{1}{D}\log \left(\dfrac{k}{m-k}\right)=\dfrac{1}{1.7}\log 4=0.815\)

２つのモデルの推定値は、ほぼ同じである。（解答参照）

(4) １母数ロジスティックモデル

識別力 \(a_j=a\) は一定、困難度 \(b_j\) が j によって異なり、

　　　　　\(p_j(\theta) =\dfrac{1}{ 1+\exp \left[ -Da\left(\theta -b_{j}\right) \right]}\)　であるとき、

θの最尤推定量を求める。尤度関数は、

\(\begin{split}L(\theta)&=\prod_{j=1}^{n} p_j(\theta)^{x_j}\left(1-p_j(\theta)\right)^{1-x_j}
=\prod_{j=1}^{n}\left(\dfrac{1}{ 1+\exp \left[ -Da\left(\theta -b_{j}\right) \right]}\right)^{x_j}
\left(\dfrac{\exp \left[ -Da\left(\theta -b_{j}\right) \right]}{ 1+\exp \left[ -Da\left(\theta -b_{j}\right) \right]}\right)^{1-x_j}\\
&=
\prod_{j=1}^{n}
\dfrac{\left(\exp \left[ -Da\left(\theta -b_{j}\right) \right]\right)^{1-x_j}}{ 1+\exp \left[ -Da\left(\theta -b_{j}\right) \right]}
\end{split}\)

対数尤度関数は、

　　　　　\(\begin{split}
l(\theta)&=\sum_{j=1}^{n}\left[-Da(\theta -b_{j})(1-x_j)\right]
-\sum_{j=1}^{n}\log\left(1+\exp \left[ -Da\left(\theta -b_{j}\right) \right]\right)
\end{split}\)

微分して　　\(\displaystyle l'(\theta)=\sum_{j=1}^{n}\left[-Da(1-x_j)\right]
-\sum_{j=1}^{n}\frac{-Da\exp \left[ -Da\left(\theta -b_{j}\right) \right]}{1+\exp \left[ -Da\left(\theta -b_{j}\right) \right]}=0\)　とすると、

　　　　　\(\displaystyle\sum_{j=1}^{n}\left(x_j-1\right)
+\sum_{j=1}^{n}\frac{\exp \left[ -Da\left(\theta -b_{j}\right) \right]}{1+\exp \left[ -Da\left(\theta -b_{j}\right) \right]}=0\)

　　　　　\(\displaystyle\sum_{j=1}^{n}x_j
=\sum_{j=1}^{n}\left[1-\frac{\exp \left[ -Da\left(\theta -b_{j}\right) \right]}{1+\exp \left[ -Da\left(\theta -b_{j}\right) \right]}\right]\)

　　　　　\(\displaystyle\therefore \sum_{j=1}^{n}x_j
=\sum_{j=1}^{n}\frac{1}{1+\exp \left[ -Da\left(\theta -b_{j}\right) \right]}\)　‥‥①

観測値 \(x_1,x_2,\cdots,x_n \) が与えられれば①を満たす \(\theta\) が決まる。

(5) ①式より、　　　\(\displaystyle\sum_{j=1}^{n}x_j
=\sum_{j=1}^{n}p_j(\theta)\)

これは、(正答数)＝(正答率の和) を意味する。

\(x_j\sim Ber(p_j(\theta))\) ベルヌーイ分布に従い、 \(E(x_j)=p_j(\theta)\) であるから、

（正答数)=(期待正答数）となる。

2015人文② 主成分分析と因子分析

簡単のため、確率変数 \(X_1,\cdots,X_p\) の期待値は 0 とする。

確率変数 \(X
=\begin{pmatrix}
X_{1}\\

\vdots\\
X_{p}
\end{pmatrix}\) の分散共分散行列 \(\displaystyle \Sigma=E(XX')\) に主成分分析を適用すると次のようになる。

（標本分散共分散行列Sに適用する場合は、解答のようになる）

新変数 \(X'w=w_1X_1+\cdots w_pX_p\) とすると、\(E(X'w)=0\) だから

分散 \(V(X'w)=E(w'XX'w)=w'E(XX')w=w'\Sigma w\)

ラグランジュの未定乗数法によると、 \(|w|=1\) のもとで \(w'\Sigma w\) を最大とするには

　　　　　\(L=w'\Sigma w-\lambda(w'w-1)\)　　を \(w\) で偏微分する。

行列と二次形式の微分を用いると、

　　　　　\(\dfrac{\partial L}{\partial w} = 2 \Sigma w - 2 \lambda w=0\)　　　　　\(\therefore \Sigma w =\lambda w\)

\(w\) は固有ベクトルであることが分かる。

\(w'\Sigma w=w'\lambda w=\lambda |w|=\lambda\) だから、固有値 \(\lambda\) は新変数の分散に等しい。

(1) \(f(w)=w'\Sigma w\) は主成分の分散であり、分散共分散行列の最大固有値に等しい。

本問では、標本分散共分散行列Sについて述べられているので、確率変数でなくn個の実データで説明する。

科目 \(X_j\) は、n人の生徒の第 \(j\) 科目の点数を成分とする列ベクトルとする。

　　　　　\(X_j=(x_{ij})
=\begin{pmatrix}
x_{1j}\\
\vdots\\
x_{ij}\\
\vdots\\
x_{nj}
\end{pmatrix}\)　　　　全ての科目の点数は \(
\begin{pmatrix}
\cdots&x_{1j}&\cdots\\
&\vdots\\
\cdots&x_{ij}&\cdots\\
&\vdots\\
\cdots&x_{nj}&\cdots
\end{pmatrix}\) （ \(n \times p\) 行列）

第 \(j\) 科目の平均を \(x_{\cdot j}\) とすると、各成分の期待値を0にでき、

　　　　　\(X=
\begin{pmatrix}
\cdots&(x_{1j}-x_{\cdot j})&\cdots\\
&\vdots\\
\cdots&(x_{ij}-x_{\cdot j})&\cdots\\
&\vdots\\
\cdots&(x_{nj}-x_{\cdot j})&\cdots
\end{pmatrix}\)　　　（列平均を０としたデータ行列）

（解説の冒頭の \(X\) とは異なるので注意）

主成分係数を \(w=(w_1,w_2,\cdots,w_p)'\) とすると、n人の主成分得点は \(Xw\)

n人の主成分得点の平方和は、 \(w'X'Xw\) だから、分散は \(f(w)=\dfrac{1}{n-1}w'X'Xw\)

標本分散共分散行列を \(S\) とすると、\(S=\dfrac{1}{n-1}X'X\) だから、　\(f(w)=w'Sw\)

ラグランジュの未定乗数法を用いると、

　　・ \(w'Sw\) を最大にする \(w\) は \(S\) の固有ベクトルであり、

　　　 \(w'Sw\) の最大値は、 \(S\) の固有値の中の最大のものであり第１主成分の分散である。

(2) 確率変数 \(X
=\begin{pmatrix}
X_{1}\\

\vdots\\
X_{p}
\end{pmatrix}\) に因子分析を行う。因子 \(F=\begin{pmatrix}
F_1\\
\vdots\\
F_m
\end{pmatrix}\) 、誤差項 \(U=\begin{pmatrix}
\varepsilon_1\\
\vdots\\
\varepsilon_m
\end{pmatrix}\) とすると、

構造方程式　 \(X=\Lambda F+U\)

　　　　　\(XX'=(\Lambda F+U)(F'\Lambda+U')
=\Lambda(FF')\Lambda+\Lambda (FU')+(UF')\Lambda +UU'\)　　　だから、

因子と誤差項は無相関と考えると、 \(E(FU')=0\:,\:E(UF')=0\)

\(UU'=\begin{pmatrix}
\varepsilon_1\\
\vdots\\
\varepsilon_m
\end{pmatrix}

\begin{pmatrix}
\varepsilon_1&\cdots&\varepsilon_m
\end{pmatrix}\) について、\(E(\varepsilon_i \varepsilon_j)=0\) 、\(E(\varepsilon_i \varepsilon_i)=e_i^2\) 、\(D=\begin{pmatrix}
e_1^2&0&0\\
0&\ddots &0\\
0&0&e_m^2
\end{pmatrix}\)

\(X\) の分散共分散行列を \(\Sigma=E(XX')\) 、因子間相関行列を \(M=E(FF')\) とすると、

　　　　　\(\Sigma=\Lambda M\Lambda'+D\)　　　　　（因子分析の基本方程式）

条件は、

　・ \(X,F\) は標準化されている

　・異なる誤差項は無相関　 \(E(\varepsilon_i \varepsilon_j)=0\:\:\:(i\neq j)\)

　・因子と誤差項は無相関　 \(E(FU')=0\:,\:E(UF')=0\)

さらに、直交モデルでは因子間相関がないから、 \(M=I_p\) （単位行列）

この場合、 \(\Sigma=\Lambda \Lambda'+D\) が因子分析の基本方程式となる。

(3) AIC（赤池情報量基準）モデルの評価に関する問題

統計Web 様　https://bellcurve.jp/statistics/glossary/380.html

AICは統計的モデルの予測性の良さを、観測値と理論値の差（残差）を用いて評価する統計量。値が小さいほど当てはまりが良い。

パラメーター \(\theta\) のもとでの尤度を \(L(x | \theta )\) 、パラメーターの個数を \(p\) とするとき、

赤池情報量基準は、

　　　　　\(AIC= -2log(L(x | \theta ))+2p=-2\left[log(L(x | \theta ))-p\right]\)　　　で与えられる。

モデルのパラメーターの個数 \(p\) が必要になる。

しかし、因子分析のモデルにおいて、構造方程式をみたす因子負荷行列 Aは一通りではないので注意が必要である。

２因子3変量の場合を例にとると、

　　　　　\(X=A\begin{pmatrix}
F \\
G
\end{pmatrix}+\varepsilon =\left\{ AQ\right\} \times \left\{ Q^{-1}\begin{pmatrix}
F \\
G
\end{pmatrix}\right\} +\varepsilon \)

\(AQ=\begin{pmatrix}
a_{x}' & b_{x}' \\
a_{y}' & b_{y}' \\
a_{z}' & b_{z}'
\end{pmatrix}\) とおくと、 \(\begin{pmatrix}
a_{x}' & a_{y}' & a_{z}' \\
b_{x}' & b_{y}' & b_{z}'
\end{pmatrix}=Q^{-1}\begin{pmatrix}
a_{x} & a_{y} & a_{z} \\
b_{x} & b_{y} & b_{z}
\end{pmatrix}\)

これは、\(Q^{-1}\) で表される一次変換によって、

　　　　　\(\left( F,G\right) \rightarrow \left( F',G'\right) \)

　　　　　\(\begin{aligned}\left( a_{x},b_{x}\right) \rightarrow \left( a_{x}',b_{x}'\right) \\
\end{aligned}\)　　　\(\left( a_{y},b_{y}\right) \rightarrow \left( a_{y}',b_{y}'\right) \)　　　\(\left( a_{z},b_{z}\right) \rightarrow \left( a_{z}',b_{z}'\right) \)

因子を新因子に移す変換 \(Q^{-1}\) によって、因子の係数ベクトルは新因子の係数ベクトルに移ることを意味する。

因子の係数ベクトルとは、変量 \(x\) のFG平面上での座標であり、

\(\begin{pmatrix}
x \\
y \\
z
\end{pmatrix}=

F
\begin{pmatrix}
a_{x}\\
a_{y}\\
a_{z}
\end{pmatrix}
+
G
\begin{pmatrix}
b_{x}\\
b_{y}\\
b_{z}
\end{pmatrix}
+
\varepsilon

\)　　　において誤差項を除外すれば、

　　　　　\(x=\begin{pmatrix}a_x&b_x\end{pmatrix}
\begin{pmatrix}
F\\G
\end{pmatrix}\) より、変量 \(x\) 　　　\(\begin{pmatrix}a_x&b_x\end{pmatrix}\)　のFG平面上での座標は

　　　　　\(y=\begin{pmatrix}a_y&b_y\end{pmatrix}
\begin{pmatrix}
F\\G
\end{pmatrix}\)

　　　　　\(z=\begin{pmatrix}a_z&b_z\end{pmatrix}
\begin{pmatrix}
F\\G
\end{pmatrix}\) より、\(z\) のFG平面上での座標は \(\begin{pmatrix}a_z&b_z\end{pmatrix}\)

というベクトルで表される。

\(Q\) が直交行列であれば、 \(Q\) は回転であり、回転により \(b_x'=0\) となるようにできる。

２因子の場合はやさしく、 \(\begin{pmatrix}a_x&b_x\end{pmatrix}\) の偏角を \(\theta\) 、 \(Q\) を \(\theta\) 回転、 \(Q^{-1}\) を \((-\theta)\) 回転とすれば、

　　　　　\(\begin{pmatrix}
a_{x}' & a_{y}' & a_{z}' \\
0 & b_{y}' & b_{z}'
\end{pmatrix}=Q^{-1}\begin{pmatrix}
a_{x} & a_{y} & a_{z} \\
b_{x} & b_{y} & b_{z}
\end{pmatrix}\)　　　　　\(\therefore AQ=\begin{pmatrix}
a_{x}' & 0 \\
a_{y}' & b_{y}' \\
a_{z}' & b_{z}'
\end{pmatrix}\)

因子負荷行列は回転の自由度を持ち、回転によって１個の成分を0に固定できるから、因子負荷行列のパラメーター個数は5、これに誤差項のパラメーター3を加えると8になる。これから、

　　　　　\(AIC= -2log(L(x | \theta ))+2p=-2\left[log(L(x | \theta ))-p\right]\)

において \(p=8\) と考えてAICを計算することができる。

さて、一般の構造方程式　\(X=\Lambda F+U\) に戻る。（この先は難しくなります）

これは、\(m\) 因子、 \(p\) 変量であるから、 \(\Sigma\) は \(p\times m\) 行列。 \(Q\) を直交行列（回転）として、

　　　　　\(X=\Lambda F+U=\Lambda(QQ') F+U=(\Lambda Q)(Q' F)+U\)

構造方程式は、直交行列（回転）だけの自由度を持つ。

２因子3変量の場合の、\(\begin{pmatrix}
a_{x}' & a_{y}' & a_{z}' \\
b_{x}' & b_{y}' & b_{z}'
\end{pmatrix}=Q^{-1}\begin{pmatrix}
a_{x} & a_{y} & a_{z} \\
b_{x} & b_{y} & b_{z}
\end{pmatrix}\) と同様に、

　　　　　因子の係数ベクトルはによって新因子の係数ベクトルに移る。

・第１成分軸が第１因子係数ベクトルに重なるように回転すれば、

　　　　　 \(Q^{-1}\) によって第１因子の係数ベクトルは \(\begin{pmatrix}
\lambda_{11},0,\cdots,0
\end{pmatrix}\) に移り、さらに

・第2成分軸が第2因子係数ベクトルに重なるように回転を合成すれば、

　　　　　 \(Q^{-1}\) によって第2因子の係数ベクトルは \(\begin{pmatrix}
\lambda_{21},\lambda_{22},0\cdots,0
\end{pmatrix}\) に移り、さらに

・同様な操作を続け、第 \(m\) 成分軸が第 \(m\) 因子係数ベクトルに重なるように回転を合成すれば、

　　　　　 \(Q^{-1}\) によって第 \(m\) 因子の係数ベクトルは \(\begin{pmatrix}
\lambda_{p1},\lambda_{p2},\cdots,\lambda_{pm}
\end{pmatrix}\) に移る。

0の個数は、 \(1+\cdots+(m-1)=\dfrac{1}{2}m(m-1)\) であるので、因子負荷行列の右上三角の \(\dfrac{1}{2}m(m-1)\) 個は0になるようにできるから、

　　　　　因子負荷行列のパラメータ数は、　 \(pm-\dfrac{1}{2}m(m-1)\)

誤差項のパラメーター \(p\) 個を加えれば、

　　　　　パラメーターの総数は　　\(p+pm-\dfrac{1}{2}m(m-1)\)

(4) ２因子5変量（因子間相関がない直交モデル）

構造方程式　\(\begin{pmatrix}
X_1 \\
X_2 \\
X_3\\
X_4\\
X_5
\end{pmatrix}=\begin{pmatrix}
\lambda_{11} & \lambda_{12} \\
\lambda_{21} & \lambda_{22} \\
\lambda_{31} & \lambda_{32} \\
\lambda_{41} & \lambda_{42} \\
\lambda_{51} & \lambda_{52}
\end{pmatrix}\begin{pmatrix}
f_1 \\
f_2
\end{pmatrix}+\begin{pmatrix}
u_1 \\
u_2\\
u_3\\
u_4 \\
u_5
\end{pmatrix}\) 　　　基本方程式　\(\Sigma=\Lambda \Lambda'+D\)

Dは、対角成分以外は0であり、対角成分には独自因子の分散が並ぶ。

\(\Sigma\) は相関行列であり、対角成分は１

　　　　　\(\Lambda \Lambda'=\begin{pmatrix}
\lambda_{11} & \lambda_{12} \\
\lambda_{21} & \lambda_{22} \\
\lambda_{31} & \lambda_{32} \\
\lambda_{41} & \lambda_{42} \\
\lambda_{51} & \lambda_{52}
\end{pmatrix}

\begin{pmatrix}
\lambda_{11} & \lambda_{21} &\lambda_{31} & \lambda_{41} &\lambda_{51} \\
\lambda_{12} & \lambda_{22} &\lambda_{32} & \lambda_{42} &\lambda_{52}
\end{pmatrix}\)　　の \((3,3)\) 成分は、\(\lambda_{31}^2+\lambda_{32}^2\) （共通性）

数学 \(X_3\) に注目すると、

　　　　　\(X_3=\lambda_{31}\cdot f_1+\lambda_{32}\cdot f_2+u_3\)

数学の分散は標準化されて１であり、因子間相関はないとしているので、両辺の分散をとると、

　　　　　\(V(X_3)=\lambda_{31}^2V(f_1)+\lambda_{32}^2V(f_2)+V(u_3)\)

　　　　　\(\therefore 1=(\lambda_{31}^2+\lambda_{32}^2)+0.48\)

0.48は誤差分散、0.52は共通性である。

これは、 \(X_3\) の52%が共通因子で説明でき、48%が誤差変動であることを意味している。

2015人文③ 多次元尺度構成法（MDS）

非類似度 \(s_{ij}\) は対象 \(i\) と対象 \(j\) の距離に相当する。

\(n\) 個の対象の非類似度行列 \(S=(s_{ij})\) は対角成分が０の \(n\times n\) 対称行列で距離行列に相当する。

この対称行列をもとに、 \(n\) 個の対象の \(p\) 次元座標を求めようとするのがMSDの目的である。

　　　　　座標行列 \(X=(x_{ik})\) とおく。 \(X\) は \(n\times p\) 行列 \((i=1,\cdots ,n\:\:,\:\:k=1,\cdots ,p)\)

\(XX'\) は \(n\times n\) 対称行列であり、成分が標準化された確率変数の場合、相関行列に相当する。

距離行列 \(S=(s_{ij})\) と相関行列 \(XX'\) をつなぎ、座標 \(X\) を求めることが本問の目的である。

\(S=XX'\) としたいところであるが、この等式は適当ではない。

そもそも、\(S\) は距離行列で対角成分が0であるが、 \(XX'\) の対角成分は0ではない。

そこで、 \(X\) から得られる距離の２乗の行列 \(D^{(2)}\) と、 \(S^{}\) の要素の２乗を成分とする \(S^{(2)}\) を比較する。

(1)(2) 距離の２乗行列 \(D^2\) の \((i,j)\) 成分を \(d_{ij}^2\) とすると、

　　　　　\(d_{ij}^2=(x_i-x_j)'(x_i-x_j)=x_i'x_i-2x_ix_j'+x_j',x_j\)

本問では、 \(X=
\begin{pmatrix}
x_1'\\
\vdots\\
x_n'
\end{pmatrix}\) \(n\times p\) 行列とおかれているので、\(XX'=\begin{pmatrix}
x_1'\\
\vdots\\
x_n'
\end{pmatrix}
\begin{pmatrix}
x_1&\cdots&x_n
\end{pmatrix}
=(x_i'x_j)\)

3次の場合で具体的に書くと、

\(diag(XX')=\begin{pmatrix}
x_1'x_1&0&0\\
0&x_2'x_2&0\\
0&0&x_3'x_3
\end{pmatrix}\)

\(diag(XX')\cdot 11'=
\begin{pmatrix}
x_1'x_1&0&0\\
0&x_2'x_2&0\\
0&0&x_3'x_3
\end{pmatrix}

\begin{pmatrix}
1&1&1\\
1&1&1\\
1&1&1
\end{pmatrix}
=
\begin{pmatrix}
x_1'x_1&x_1'x_1&x_1'x_1\\
x_2'x_2&x_2'x_2&x_2'x_2\\
x_3'x_3&x_3'x_3&x_3'x_3
\end{pmatrix}\)

\(11'\cdot diag(XX')=
\begin{pmatrix}
1&1&1\\
1&1&1\\
1&1&1
\end{pmatrix}
\begin{pmatrix}
x_1'x_1&0&0\\
0&x_2'x_2&0\\
0&0&x_3'x_3
\end{pmatrix}
=
\begin{pmatrix}
x_1'x_1&x_2'x_2&x_3'x_3\\
x_1'x_1&x_2'x_2&x_3'x_3\\
x_1'x_1&x_2'x_2&x_3'x_3
\end{pmatrix}\)

\(XX'=\begin{pmatrix}
x_1'\\
x_2'\\
x_3'
\end{pmatrix}
\begin{pmatrix}
x_1&x_2&x_3
\end{pmatrix}
=
\begin{pmatrix}
x_1'x_1&x_1'x_2&x_1'x_3\\
x_2'x_1&x_2'x_2&x_1'x_3\\
x_3'x_1&x_3'x_2&x_1'x_3
\end{pmatrix}\)

よって、 \(D^{(2)}=diag(XX')\cdot 11'-2XX'+11'\cdot diag(XX')\)

\(S^{(2)}=D^{(2)}\) となる \(X\) を求めたいので、２重中心化したものを比較する。

［２重中心化行列］

中心化行列は、\(J_n=I_n-\dfrac{1}{n}1_n\cdot 1_n'\) （ \(J_n\) を左からかければ列平均は0、右からかければ行平均は0）

各成分の平均を0にするには、\(X=(x_{ik}-x_{\cdot k})\) を考えればよいが、

全ての成分を1とするベクトルを１とすると、 \(1\cdot 1'=
\begin{pmatrix}
1&1&\cdots & 1\\
1&1&\cdots & 1\\
\vdots&\vdots& \ddots & \vdots\\
1&1&\cdots & 1
\end{pmatrix}
\) だから

各成分の平均は、\(\dfrac{1}{n}1\cdot 1'X\) で表されるので、列平均を0とするには

中心化行列を \(J_n=I_n-\dfrac{1}{n}1\cdot 1'X\) とおいて、 \(J_nX\) を計算すればよい。

\(X\) の列平均が０である条件は、 \(X=J_nX\) である。

次に、二重中心化とは、中心化行列を左右からかけて \(-\dfrac{1}{2}\) 倍する操作を言う。

　　　　　\(-\dfrac{1}{2}J_n\cdot S^{(2)}\cdot J_n=-\dfrac{1}{2}J_n\cdot D^{(2)}\cdot J_n\)　　　あるいは \(-\dfrac{1}{2}J_n\cdot S^{(2)}\cdot J_n \fallingdotseq -\dfrac{1}{2}J_n\cdot D^{(2)}\cdot J_n\)

を満たす \(X\) を求めることを目的としよう。

ここで、\(1'J_n=0'\) 、\(J_n1=0\) だから、\(J_nD^{(2)}J_n=-2J_nXX'J_n=-2(J_nX)(J_nX)'\)

\(X\) を中心化（各成分の平均を0）した座標を \(X^*=J_nX\) とすると、 \(-\dfrac{1}{2}J_nD^{(2)}J_n=X^{*} (X^{*})'\)

\(X\) はあらかじめ中心化してあるとすると、

\(X=J_nX=X^*\) であるので、 \(-\dfrac{1}{2}J_nD^{(2)}J_n=XX'\)

よって、　\(-\dfrac{1}{2}J\cdot S^{(2)}\cdot J=XX'\) 　を満たす \(X\) を求める問題に帰着する。

(3) \(B=-\dfrac{1}{2}J\cdot S^{(2)}\cdot J\) とおいて、 \(B=XX'\) を求める問題を考える。

中心化された対称行列 \(B\) に対し \(B=XX'\) を満たす \(X\) を求める問題は、因子分析で因子負荷行列を求める問題と同様の手法が利用できる。

\(B\) は \(n\times n\) 対称行列であるから、直交行列 \(Q\) で対角化でき、

　　　　　\(Q^{-1} BQ=\Lambda\)　　　　　\(\therefore B=Q\Lambda Q^{-1}\)

　　　　　\(XX'=Q\Lambda Q^{-1}\)

もし、固有値が0以上ならば、固有値の正の平方根を対角成分とする対角行列を \(\Lambda^{1/2}\) とすれば、

　　　　　\(XX'=\left(Q\Lambda^{1/2}\right)\left(\Lambda^{1/2} Q^{-1}\right)=\left(Q\Lambda^{1/2}\right)\left(Q\Lambda^{1/2}\right)'\)

\(X=Q \Lambda^{1/2}\) としたいところであるが、\(Q\Lambda^{1/2}\) は \(n\times n\) 行列、 \(X\) は \(n\times p\) 行列である。

\(B\) のスペクトル分解　\(\displaystyle B=\sum_{i=1}^n \lambda_iq_iq_i'\) で絶対値の小さい固有値と固有ベクトルを無視しても、

　　　　　\(\displaystyle B\fallingdotseq \sum_{i=1}^p \lambda_iq_iq_i\)　　　が近似的に言える。

成分で具体的に説明する。固有値を降順に対角に配置しておく。

\(X\) は \(n\times p\) 行列である。上の \(Q\) は \(n\times n\) 行列であるが、はじめの \(p\) 列 \(Q_+\) は、 \(n\times p\) 行列である。

\(\Lambda\) のはじめの \(p\times p \) 小行列 \(\Lambda_+\) は、 \(p\times p\) 行列であり、

　　　　　\(Q_+^{-1} BQ_+=\Lambda_+\)　　　　　次数は \((p,n)\times (n,n)\times (n,p)=(p,p)\)

　　　　　\(XX'=Q_+\Lambda_+ Q_+^{-1}\)　　　　次数は \((n,p)\times (p,n)=(n,p)\times (p,p)\times (p,n)\)

\(\Lambda_+=\begin{pmatrix}
\lambda_{1}&0&0\\
0&\ddots &0\\
0&0&\lambda_p
\end{pmatrix}\) とおくと、\(\Lambda_+^{1/2}=\begin{pmatrix}
\sqrt{\lambda_{1}}&0&0\\
0&\ddots &0\\
0&0&\sqrt{\lambda_p}
\end{pmatrix}\) 、 \(\Lambda_+=\Lambda_+^{1/2}\cdot \Lambda_+^{1/2}\) だから、

　　　　　\(XX'\fallingdotseq \left(Q_+\Lambda_+^{1/2}\right)\left(\Lambda_+^{1/2} Q_+^{-1}\right)\)　　であり、 \(X\) も \(Q_+\Lambda_+^{1/2}\) も \(n\times p\) 行列だから、

\(X\fallingdotseq Q_+\Lambda_+^{1/2}\) とすればよいことが分かる。

(4) \(U\) を直交行列とすると、 \(XX'=XUU'X'\) が成立するから、 \(XU\) も解である。

つまり因子分析の場合と同様に、\(X\) は回転だけの自由度を持つ。

近似について、もっと詳しく考えてみよう。

(3)について、実際には、 \(X\) は近似解で求まればよいので、

　　　　　\(\left \|-\dfrac{1}{2}J_n\cdot S^{(2)}\cdot J_n -\left( -\dfrac{1}{2}J_n\cdot D^{(2)}\cdot J_n\right)\right \|\)　　　を最小とする \(X\) を求めることになる。

そもそも、MSDは次元圧縮を伴うので、数学的な唯一解ではなく近似解である。

　　　　　\(f=\left \|B-XX'\right \|^2\)　　　を最小とする \(X\) を微分を利用して求める。

行列の微分では、 \(\dfrac{\partial}{\partial X}XX'=2X\) 　　だから、　\(\dfrac{\partial}{\partial x}f=-4(B-XX')X\)

\(\dfrac{\partial}{\partial x}f=0\) とすると、\(BX-X(X'X)=0\)

MSDの目的は、距離行列から座標行列 \(X\) を算出することにより距離情報を可視化することであるから、 \(X\) には回転だけの自由度がある。

\(X'X\) は対称行列であるから、 \(P'(X'X)P\) が対角行列になるように、直交行列 \(P\) を選び \(XP\) をあらためて \(X\) と考えれば、 \(X'X=\Delta\) が対角行列になるようにできる。

　　　　　\(BX-X \Delta =0\)　より、\(BX=X \Delta\)

\(X\) の列ベクトルが \(B\) の固有ベクトル、\(\Delta\) の対角成分が固有値となる。

［フロベニウスノルム］　　　 \(\|A\|_F^2 = Tr(AA^\top)\) を利用する。

たふぇたふぇのブログ様　参照
http://blog.tafetafe.work/studies/mathematics/linear-algebra/frobenius-norm-and-trace/#outline__2

トレースは線型性、可換性を持つことに注意する。

　　　　　\(\begin{split}
f&=\left \|B-XX'\right \|^2=Tr\left((B-XX')(B'-XX') \right)\\
&=Tr\left(BB'\right)+Tr\left(XX'XX'\right)-Tr\left(BXX' \right)-Tr\left(XX'B' \right)\\
&=Tr\left(B^2\right)+Tr\left(XX'XX'\right)-2Tr\left(BXX' \right)\\
&=Tr\left(B^2\right)+Tr\left(\Delta^2\right)-2Tr\left(X \Delta X' \right)\\
&=Tr\left(B^2\right)+Tr\left(\Delta^2\right)-2Tr\left( \Delta X'X\right)\\
&=Tr\left(B^2\right)+Tr\left(\Delta^2\right)-2Tr\left(\Delta^2\right)\\
&=Tr\left(B^2\right)-Tr\left(\Delta^2\right)
\end{split}\)

\(B\) の固有値の平方和が最大のとき、 \(f\) は最小となることが分かる。

つまり、固有値を絶対値の大きい方から \(p\) 個選んで、 \(X\fallingdotseq Q_+\Lambda_+^{1/2}\) とすれば、

\(B\) と\(XX'\) のフロベニウスのノルム（行列の距離）が最小となることが分かる。

［参考］

初級Mathマニアの寝言様
　　　　　https://ogyahogya.hatenablog.com/entry/2019/12/30/特異値分解

2015人文④ ベイズ予測ルール

ベイズ推定は、得られている情報から未知の状況を推定するものである。

財務指標 \(X\) を説明変数として、倒産するかどうかの反応変数 \(Y\) を予測する。

（ \(X\) は情報、 \(Y=1,2\) はそれぞれ倒産、非倒産という状況を表す）

倒産した会社について、１０年前の財務指標 \(X\) を調査すれば、 \(X\) の統計情報を得ることは可能である。

本問では、 \(Y=1\) のとき \(X\sim N(10,10^2)\) 、 \(Y=2\) のとき \(X\sim N(40,10^2)\) が分かっているものとする。

事前確率は、 \(P(Y=1)=\dfrac{1}{2}\:,\:P(Y=2)=\dfrac{1}{2}\)

\(x\) の尤度は、 \(P(X=x|H=1)\:,\,P(X=x|H=2)\) であり、

　　　　　\(P(X=x|H=1)=\dfrac{1}{\sqrt{2\pi }\cdot 10}\exp \left( -\dfrac{(x-10)^{2}}{2}\right)\)

　　　　　\(P(X=x|H=2)=\dfrac{1}{\sqrt{2\pi }\cdot 10}\exp \left( -\dfrac{(x-40)^{2}}{2}\right)\) 　　　　で与えられる。

\(P(Y=1|X=x)\propto P(H=1)\cdot P(X=x|H=1)\) 、\(P(Y=2|X=x)\propto P(H=2)\cdot P(X=x|H=2)\)

事後確率は条件付き確率だから、確率にするには正規化して合計を１にする必要がある。

しかし、事後確率 \({P(Y=2|X=x)}\:\:,\:\:{P(Y=1|X=x)}\) の大小比較だけなら、割り算して１と比較すればよい。

　　　　　\(\dfrac{P(Y=2|X=x)}{P(Y=1|X=x)}=\dfrac{P(H=2)\cdot P(X=x|H=2)}{P(H=1)\cdot P(X=x|H=1)}=\exp \left( -\dfrac{(x-40)^{2}}{2}+\dfrac{(x-10)^{2}}{2}\right)\)

これが１より大きいかどうか判断するには、対数をとって正か負かを調べればよい。

　　　　　\(-\dfrac{(x-40)^{2}}{2}+\dfrac{(x-10)^{2}}{2}=\dfrac{(2x-50)(30)}{2}=30(x-25)\)

よって、 \(x-25>0\) のとき　\(P(Y=2|X=x)>P(Y=1|X=x)\)

倒産していない確率の方が大きいと判断することができる。

(2) 倒産しているのに倒産していないと判断する確率は、

　　　　　\(\displaystyle P(X>25|H=1)=\int _{25}^{+\infty}\dfrac{1}{\sqrt{2\pi }\cdot 10}\exp \left( -\dfrac{(x-10)^{2}}{2}\right)dx\)

\(X=25\) は標準化すると、\(Z=1.5\)

正規分布表より、 \(\displaystyle \int _{1.5}^{+\infty}\dfrac{1}{\sqrt{2\pi }}\exp \left( -\dfrac{z^{2}}{2}\right)dz=0.0668\) 　　　　 \(\therefore P(X>25|H=1)=0.0068\)

(3) 説明変数である財務指標が、\(X=\begin{pmatrix}
X_1\\
X_2
\end{pmatrix}\) となった場合を考える。

\(X_1,X_2\) は独立ではなく、2変量正規分布に従っていると考えるのが自然である。

　　　　　\(Y=1\)　のとき \(X=\begin{pmatrix}
X_1\\
X_2
\end{pmatrix}
\sim N\left(\mu^{(1)},\Sigma\right)\) 、 \(Y=2\) のとき \(X=\begin{pmatrix}
X_1\\
X_2
\end{pmatrix}
\sim N\left(\mu^{(2)},\Sigma\right)\)

分散共分散行列は共通であるとする。

\(x=\begin{pmatrix}
x_1\\
x_2
\end{pmatrix}\) の尤度は、 \(P(X=x|H=1)=\dfrac{1}{\left(\sqrt{2\pi}\right)^2\sqrt{|\Sigma|}}\exp\left[{-\dfrac{1}{2}(x-\mu^{(1)})'\ \Sigma^{-1}(x-\mu^{(1)})}\right]\)

　　　　　　　　　\(P(X=x|H=2)=\dfrac{1}{\left(\sqrt{2\pi}\right)^2\sqrt{|\Sigma|}}\exp\left[{-\dfrac{1}{2}(x-\mu^{(2)})'\ \Sigma^{-1}(x-\mu^{(2)})}\right]\)

　　　　　\(\therefore \dfrac{P(Y=2|X=x)}{P(Y=1|X=x)}=\exp\left[{-\dfrac{1}{2}(x-\mu^{(2)})'\ \Sigma^{-1}(x-\mu^{(2)})}+\dfrac{1}{2}(x-\mu^{(1)})'\ \Sigma^{-1}(x-\mu^{(1)})\right]\)

対数をとると、

　　　　　\(\begin{split}
-\dfrac{1}{2}(x-\mu^{(2)})'\ \Sigma^{-1}(x-\mu^{(2)})+\dfrac{1}{2}(x-\mu^{(1)})'\ \Sigma^{-1}(x-\mu^{(1)})

=
\left(\mu^{(2)}-\mu^{(1)}\right)'\ \Sigma^{-1}\left(x-\dfrac{\mu^{(1)}+\mu^{(2)}}{2}\right)\end{split}\)

\(\left(\mu^{(2)}-\mu^{(1)}\right)'\ \Sigma^{-1}\left(x-\dfrac{\mu^{(1)}+\mu^{(2)}}{2}\right)=0\) は、\(\mu^{(1)},\mu^{(2)}\) の中点を通り \(\left(\mu^{(2)}-\mu^{(1)}\right)'\ \Sigma^{-1}\) を法線ベクトルとする直線、即ち、線形関数である。

（これは線形判別式を求める計算と同じである。当サイト「線形判別」参照）

(4) ベイズ予測ルールを題材としているが、判別分析と同様の計算で判別境界を求めている。

\(X_1,X_2\) の2変量正規分布の分散共分散行列は、相関係数を \(\rho \) とすると、

　　　　　\(\Sigma=\begin{pmatrix}
\sigma _{1}^{2} & \rho\sigma_{1}\sigma _{2} \\
\rho\sigma_{1}\sigma _{2} & \sigma _{2}^{2}
\end{pmatrix}\)　　　　　\(\therefore \Sigma^{-1}=\dfrac{1}{\sigma_1^2\sigma_2^2(1-\rho^2)}
\begin{pmatrix}
\sigma _{2}^{2} & -\rho\sigma_{1}\sigma _{2} \\
-\rho\sigma_{1}\sigma _{2} & \sigma _{1}^{2}
\end{pmatrix}\)

\(\sigma_1=\sigma_2\) だから、\(\Sigma^{-1}=\dfrac{1}{\sigma^2(1-\rho^2)}
\begin{pmatrix}
1 & -\rho \\
-\rho &1 \end{pmatrix}
\)

　　　　　\(\left(\mu^{(2)}-\mu^{(1)}\right)'\ \Sigma^{-1}=\begin{pmatrix}30&30\end{pmatrix}
\dfrac{1}{10^2(1-\rho^2)}
\begin{pmatrix}
1 & -\rho \\
-\rho & 1
\end{pmatrix}
=\dfrac{30}{10^2(1-\rho^2)}
\begin{pmatrix}
1-\rho & 1-\rho
\end{pmatrix}\)

\(\mu^{(1)},\mu^{(2)}\) の中点 \((25,25)\) を通り \((1,1)\) を法線ベクトルとする直線は、

　　　　　\((1,1)\cdot (x_1-25,x_2-25)=0\)　　　　　\(\therefore x_1+x_2=50\)

倒産していないと判断されるのは、 \((1,1)\cdot (x_1-25,x_2-25)>0\) のときであるから、

倒産していて倒産していないと判断される確率は、 \(P(X_1+X_2>50|Y=1)\) である。

\(Y=1\) のとき \(X=
\begin{pmatrix}
X_1 \\
Y_2
\end{pmatrix}\sim N\left( \begin{pmatrix}
10 \\
10
\end{pmatrix},
\begin{pmatrix}
10^2 & 10^2\rho\\
10^2\rho & 10^2
\end{pmatrix}\right)
\) だから、

\(X_1+X_2\) の分散は、

　　　　　\(\begin{split}
V(X_1+X_2)&=V(x_1)+V(X_2)+2cov(X_1,X_2)\\
&=V(x_1)+V(X_2)+2\dfrac{cov(X_1,X_2)}{\sqrt{V(X_1)}\sqrt{V(X_2)}}\sqrt{V(X_1)}\sqrt{V(X_2)}\\
&=V(x_1)+V(X_2)+2\rho\sqrt{V(X_1)}\sqrt{V(X_2)}
\end{split}\)

\(X_1+X_2 \sim N\left(20\:,\:2\sigma^2(1+\rho)\right)= N\left(20\:,\:2\cdot10^2\dfrac{9}{8}\right)=N(20,15^2)\) だから、

\(\displaystyle P\left(X_1+X_2>50|Y=1\right)=P\left(\dfrac{X_1+X_2-20}{15}>2|Y=1\right)=\int _{2}^{+\infty}\dfrac{1}{\sqrt{2\pi }}\exp \left( -\dfrac{z^{2}}{2}\right)dz=0.0228\)

2015人文⑤ ２元配置分散分析

東京大学大学院大森宏様　「分散分析」
http://lbm.ab.a.u-tokyo.ac.jp/~omori/meiji/meiji04new.html
を参考にさせていただきました。

(1)(2)(3)は解答参照。(4)について説明する。

偏差の性質　 \(\displaystyle \sum_n(x_n-\overline{x})=0\) （平均を引いておけば、合計は0）

一般の \(a\) に対し、 \(\displaystyle \sum_n(x_n-a)^2
=n(\overline{x}-a)^2
+\sum_n(x_n-\overline{x})^2\) 　が成立する。

\((A_i,B_j)\) 水準に10個のデータ \(x_{ijk}\:\:(k=1:10)\) があるとすると、

　　　　　\(\displaystyle \sum_{k=1}^{10}(x_{ijk}-a)^2
=10(\overline{x}_{ij\cdot}-a)^2
+\sum_{k=1}^{10}(x_{ijk}-\overline{x}_{ij\cdot})^2\)

\(a=\overline{x}_{\cdot\cdot\cdot}\) とすると、

　　　　　\(\displaystyle \sum_{k=1}^{10}(x_{ijk}-\overline{x}_{\cdot\cdot\cdot})^2
=10(\overline{x}_{ij\cdot}-\overline{x}_{\cdot\cdot\cdot})^2
+\sum_{k=1}^{10}(x_{ijk}-\overline{x}_{ij\cdot})^2\)　

これは、(偏差平方和)=(群間変動)+(群内変動) を意味する。

さらに \(i,j\) について和をとると、

　　　　　\(\displaystyle \sum_{i=1}^2\sum_{j=1}^3\sum_{k=1}^{10}(x_{ijk}-\overline{x}_{\cdot\cdot\cdot})^2
=10\sum_{i=1}^2\sum_{j=1}^3(\overline{x}_{ij\cdot}-\overline{x}_{\cdot\cdot\cdot})^2
+\sum_{i=1}^2\sum_{j=1}^3\sum_{k=1}^{10}(x_{ijk}-\overline{x}_{ij\cdot})^2\)

\(\displaystyle \sum_{i=1}^2\sum_{j=1}^3(\overline{x}_{ij\cdot}-\overline{x}_{\cdot\cdot\cdot})^2
\) をさらに平方和に分解することを考える。

　　　　　 \(A_i\) 水準の平均 \(\overline{x}_{i\cdot\cdot}\) と全平均 \(\overline{x}_{\cdot\cdot\cdot}\) の差は、\(\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot\cdot\cdot}\)

　　　　　 \(B_j\) 水準の平均 \(\overline{x}_{\cdot j\cdot}\) と全平均 \(\overline{x}_{\cdot\cdot\cdot}\) の差は、\(\overline{x}_{\cdot j\cdot}-\overline{x}_{\cdot\cdot\cdot}\)

３平方の定理のように、　\(\displaystyle \sum_{i=1}^2\sum_{j=1}^3(\overline{x}_{ij\cdot}-\overline{x}_{\cdot\cdot\cdot})^2
=\sum_{i=1}^2\sum_{j=1}^3(\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot\cdot\cdot})^2+\sum_{i=1}^2\sum_{j=1}^3(\overline{x}_{\cdot j\cdot}-\overline{x}_{\cdot\cdot\cdot})^2
\)

となってほしいところであるが、これは成立しない。

　　　　　\((\overline{x}_{ij\cdot}-\overline{x}_{\cdot \cdot\cdot})
=(\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot \cdot\cdot})
+(\overline{x}_{\cdot j\cdot}-\overline{x}_{\cdot\cdot\cdot})
+(\overline{x}_{ij\cdot}-\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot j\cdot}+\overline{x}_{\cdot \cdot\cdot})\)

　　　　　(全効果)= (Aの効果)+ (Bの効果)+ (ABの交互作用)

このように、ABの交互作用の項が入ってくる。

一般の場合で説明する。偏差の合計は0になるので、

　　　　　\(\displaystyle \sum_i(\overline{x}_{ij\cdot}-\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot j\cdot}+\overline{x}_{\cdot \cdot\cdot})=\sum_i(\overline{x}_{ij\cdot}-\overline{x}_{\cdot j\cdot})-\sum_i(\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot \cdot\cdot})=0\)

　　　　　\(\displaystyle \sum_j(\overline{x}_{ij\cdot}-\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot j\cdot}+\overline{x}_{\cdot \cdot\cdot})
=\sum_j(\overline{x}_{ij\cdot}-\overline{x}_{i\cdot\cdot})-\sum_j(\overline{x}_{\cdot j\cdot}-\overline{x}_{\cdot \cdot\cdot})=0\)

　　　　　\(\displaystyle \sum_{i}(\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot \cdot\cdot})=0\) 、 \(\displaystyle \sum_{j}(\overline{x}_{\cdot j\cdot}-\overline{x}_{\cdot\cdot\cdot})=0\) 　であるので、

\((\overline{x}_{ij\cdot}-\overline{x}_{\cdot \cdot\cdot})^2
=\left[(\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot \cdot\cdot})
+(\overline{x}_{\cdot j\cdot}-\overline{x}_{\cdot\cdot\cdot})
+(\overline{x}_{ij\cdot}-\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot j\cdot}+\overline{x}_{\cdot \cdot\cdot})\right]^2\) を展開して合計した項において、

　　　　　\(\displaystyle \sum_{ij}(\overline{x}_{\cdot j\cdot}-\overline{x}_{\cdot\cdot\cdot})
(\overline{x}_{ij\cdot}-\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot j\cdot}+\overline{x}_{\cdot \cdot\cdot})
=\sum_{j}\left[(\overline{x}_{\cdot j\cdot}-\overline{x}_{\cdot\cdot\cdot})
\sum_i(\overline{x}_{ij\cdot}-\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot j\cdot}+\overline{x}_{\cdot \cdot\cdot})\right]
=0\)

　　　　　\(\displaystyle \sum_{ij}(\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot \cdot\cdot})
(\overline{x}_{ij\cdot}-\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot j\cdot}+\overline{x}_{\cdot \cdot\cdot})
=\sum_{i}\left[(\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot \cdot\cdot})
\sum_j(\overline{x}_{ij\cdot}-\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot j\cdot}+\overline{x}_{\cdot \cdot\cdot})\right]
=0\)

　　　　　\(\displaystyle \sum_{ij}(\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot \cdot\cdot})
(\overline{x}_{\cdot j\cdot}-\overline{x}_{\cdot\cdot\cdot})
=
\sum_{i}(\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot \cdot\cdot})
\sum_{j}(\overline{x}_{\cdot j\cdot}-\overline{x}_{\cdot\cdot\cdot})
=0\)

よって、

　　　　　\(\displaystyle \begin{split}
\sum_{ij}(\overline{x}_{ij\cdot}-\overline{x}_{\cdot \cdot\cdot})^2
&=\sum_{ij}(\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot \cdot\cdot})^2
+\sum_{ij}(\overline{x}_{\cdot j\cdot}-\overline{x}_{\cdot\cdot\cdot})^2
+\sum_{ij}(\overline{x}_{ij\cdot}-\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot j\cdot}+\overline{x}_{\cdot \cdot\cdot})^2\\
\end{split}\)

即ち、 \(\displaystyle \sum_{i=1}^2\sum_{j=1}^3(\overline{x}_{ij\cdot}-\overline{x}_{\cdot\cdot\cdot})^2
=3\sum_{i=1}^2(\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot \cdot\cdot})^2
+2\sum_{j=1}^3(\overline{x}_{\cdot j\cdot}-\overline{x}_{\cdot\cdot\cdot})^2
+\sum_{i=1}^2\sum_{j=1}^3(\overline{x}_{ij\cdot}-\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot j\cdot}+\overline{x}_{\cdot \cdot\cdot})^2
\)

\(\displaystyle \begin{split}
\sum_{i=1}^2\sum_{j=1}^3\sum_{k=1}^{10}(x_{ijk}-\overline{x}_{\cdot\cdot\cdot})^2

=&30\sum_{i=1}^2(\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot\cdot\cdot})^2+20\sum_{j=1}^3(\overline{x}_{\cdot j\cdot}-\overline{x}_{\cdot\cdot\cdot})^2+
10\sum_{i=1}^2\sum_{j=1}^3(\overline{x}_{ij\cdot}-\overline{x}_{i\cdot\cdot}-\overline{x}_{\cdot j\cdot}+\overline{x}_{\cdot\cdot\cdot})^2\\
&+\sum_{i=1}^2\sum_{j=1}^3\sum_{k=1}^{10}(x_{ijk}-\overline{x}_{ij\cdot})^2
\end{split}\)

(5) 推定量の基本を確認（次のサイトを参考にさせていただいた）

熊本大学数理科学総合教育センター様　「不偏推定量, 一致推定量, 最尤推定量」
　　　　　　http://msec.kumamoto-u.ac.jp/problem/pdf/statistics/2_3/ans/ex_s2_3_1_ans.pdf
東洋大学経済学部三平剛様　「最小二乗法」
　　　　　　http://www2.toyo.ac.jp/~mihira/keizaitoukei2015/ols2.pdf
名大の授業様
　　　　　　https://ocw.nagoya-u.jp/files/247/05_mle.pdf
Shingoの数学ノート様　「最尤推定量の統計的性質」
　　　　　　http://mathshingo.chillout.jp/blog23.html

不偏推定量
　・最小二乗法の損失関数が最小となる推定値。
　・1回1回の推定値は真の値から誤差が生じるが、誤差の生じ方に偏りがなく推定される。
　・推定量の期待値が真のパラメター値に等しい。
　・撹乱項の分散の不偏推定量は、RSS を自由度で割って得られる。

最尤推定量
　・尤度関数が最大となる推定値。
　・最尤推定量 θ􏰇の分布は真のパラメータθを期待値とする正規分布に法則収束。

一致推定量
　・データ(標本)の数が増えると、推定値は真の値に限りなく収束する。
　・適当な条件のもとで最尤推定量は一致推定量である。

効率的推定量
　・推定値のバラツキ(分散)が小さく、高い精度で推定できる。（最小分散不偏推定量など）
　・データの分布により判断する必要がある。（三平先生のp7を参照）

［クラメール・ラオの下限］

不偏推定量の場合、推定量の分散はフィッシャー情報量の逆数より小さくなることはない。

　　　　　\(V[\hat{\theta}]\geq J_n(\theta)^{-1}\) 　　　　　 \(\hat{\theta}\) は不偏推定量、 \(J_n(\theta)\) はフィッシャー情報量

（参考サイト）Avilen AI Trend 様　https://ai-trend.jp/basic-study/estimator/cramer-rao/

最尤推定量は一般にパラメータの不偏推定量ではないが、

　n → ∞ という漸近的な状況で、推定量の漸近分散がフィッシャー情報量の逆数に一致する。

（参考サイト）Goodな生活様「最尤推定量とクラメール・ラオの不等式」
　　　　　　　　https://www.goodnalife.com/entry/2020/11/22/000000

計算については、解答を参照。

hyoccori-toukei

サイトナビゲーション[スキップ]

2015 人文　（統計検定1級応用　過去問の解説です）

2015人文① 項目反応理論

2015人文② 主成分分析と因子分析

2015人文③ 多次元尺度構成法（MDS）

2015人文④ ベイズ予測ルール

2015人文⑤ ２元配置分散分析

hyoccori-toukei

サイトナビゲーション[スキップ]

2015 人文 （統計検定1級応用 過去問の解説です）

2015人文① 項目反応理論

2015人文② 主成分分析と因子分析

2015人文③ 多次元尺度構成法（MDS）

2015人文④ ベイズ予測ルール

2015人文⑤ ２元配置分散分析

2015 人文　（統計検定1級応用　過去問の解説です）