1 問題1_混合分布.

1.1 (1).

偏差値の問題は解答参照.

合格率  \(P(G_1)\cdot P(X>70|G_1)+P(G_2)\cdot P(X>70|G_2)\).

  \(P(X>70|G_1)=P\left(\dfrac{X-70}{4}>0~|~G_1\right)=P(Z>0)=0.5\).
  \(P(X>70|G_2)=P\left(\dfrac{X-79}{4}>\dfrac{9}{4}~|~G_2\right)=P(Z>2.5)=1-0.0062=0.9938\).
よって、
合格率: \(\dfrac{2}{3}\cdot 0.5+\dfrac{1}{3}\cdot 0.9938=0.6646\).

1.2 (2).

平均は加重平均で計算する。
   \(\dfrac{2}{3}\cdot 0.5+\dfrac{1}{3}\cdot 0.9938=0.6646\).

全変動=群内変動+群間変動 により.

   全変動=\((200\cdot 4^2+100\cdot 4^2)+\left[200\cdot (70-73)^2+100\cdot (79-73)^2\right]\).

300で割って.

   分散=\(\left[\dfrac{2}{3}\cdot 4^2+\dfrac{1}{3}\cdot 4^2\right]+\left[\dfrac{2}{3}\cdot (70-73)^2+\dfrac{1}{3}\cdot (79-73)^2\right]=16+\dfrac{54}{3}=34\)

偏差値の問題は解答参照.

1.3 (3).

文系 \(G_1\) の確率密度関数を \(f_1(x)=\dfrac{1}{\sqrt{2\pi}\cdot 4}\exp\left[-\dfrac{(x-70)^2}{2\cdot 4^2}\right]\).

理系 \(G_2\) の確率密度関数を \(f_2(x)=\dfrac{1}{\sqrt{2\pi}\cdot 4}\exp\left[-\dfrac{(x-79)^2}{2\cdot 4^2}\right]\).

混合比を \(p_1=P(G_1)\:,\:p_2=P(G_2)\) とすると、混合分布の確率密度関数は.

   \(f(x)=p_1f_1(x)+p_2f_2(x)=\dfrac{2}{3}f_1(x)+\dfrac{1}{3}f_2(x)\).
    (このままでは確率密度関数にならない)
全積分は.
 \(\displaystyle\int f(x)dx=\dfrac{2}{3}\int f_1(x)dx+\dfrac{1}{3}\int f_2(x)dx=\dfrac{2}{3}\sqrt{2\pi}\cdot 4+\dfrac{1}{3}\sqrt{2\pi}\cdot 4=4\sqrt{2\pi}\)  だから、

   \(f(x)=\dfrac{1}{4\sqrt{2\pi}}\left[\dfrac{2}{3}f_1(x)+\dfrac{1}{3}f_2(x)\right]\) とすると 全積分は1になるので確率密度関数と言える。

1.4 (3-1).

手書きでグラフを描くのには戸惑った。
和のグラフが2峰か1峰かという問題があるが、グラフを継ぎ足して加算するしかないと思う。
グラフを継ぎ足して描けば、

  • 1峰であること.
  • \(G_1\)の極大となるxが、和の分布の極大点.

であることが分かり、微分して確かめるのは時間の無駄と思われる。

f=function (m,s,x) (1/(sqrt(2*pi)*s)*exp(-(x-m)^2/(2*(s^2))))
par(mfrow=c(1,2))
curve(2*f(70,4,x),xlim=c(65,85))
curve(f(79,4,x),xlim=c(65,85),add=T)
curve(2*f(70,4,x)+f(79,4,x),xlim=c(65,85))

1.5 (3-2).

この問題は過去に類題があり、易しいと考えて微分計算を付け足したのだと思う。
参考までに、Rのグラフで導関数の符号を調べてみると次のグラフのようになる。

curve(2*(-2)*(x-70)/32*exp(-((x-70)^2)/32)+(-2)*(x-79)/32*exp(-((x-79)^2)/32),xlim=c(65,85))
abline(h=0,col="lightblue")

2 問題2_偏相関係数とパス解析.

偏相関係数と偏回帰係数.
小問[1]は代入し計算するだけである。
小問[2]は、以下の定義と基本事項に関する問題である。

2.1 偏相関係数の定義.

偏相関係数とは、2変数の相関から他の変数の影響を除いたものである。
具体的には、
  2変数の各々を他の変数で回帰したときの残差同士の相関ということになる。

本問は3変数の場合であり、
  X,Y から共変量Wの影響を除いた偏相関係数を考える。

<XをWで回帰>.
単回帰であり、変数は標準化されているので、回帰係数=相関係数 である。
   \(X=r_{XW}\cdot W+E_1\)   残差:\(E_1=X-r_{XW}\cdot W\).

<YをWで回帰>.
   \(Y=r_{YW}\cdot W+E_2\)   残差:\(E_2=Y-r_{YW}\cdot W\).

偏相関係数は \(Cor(E_1,E_2)=\dfrac{Cov(E_1,E_2)}{\sqrt{V(E_1)}\sqrt{V(E_2)}}\).

  \(\begin{split} V(E_1)&=V(X-r_{XW}\cdot W)=V(X)-r_{XW}^2V(W)=1-r_{XW}^2 \end{split}\).
  \(\begin{split} V(E_2)=1-r_{YW}^2 \end{split}\).
  \(\begin{split} Cov(E_1,E_2)&=E(E_1E_2)=E(XY-r_{XW}\cdot YW-r_{YW}\cdot XW+r_{XW}r_{YW}\cdot W^2)\\ &=r_{XY}-r_{XW}r_{YW}-r_{YW}r_{XW}+r_{XW}r_{YW}=r_{XY}-r_{XW}r_{YW} \end{split}\).
よって、
  \(Cor(E_1,E_2)=\dfrac{r_{XY}-r_{XW}r_{YW}}{\sqrt{1-r_{XW}^2}\sqrt{1-r_{YW}^2}}\)   \(\therefore r_{XY/W}=\dfrac{r_{XY}-r_{XW}r_{YW}}{\sqrt{1-r_{XW}^2}\sqrt{1-r_{YW}^2}}\).

<確認>.
3変量において、

  • 偏相関係数は相関係数から計算される.
  • 偏相関係数は、回帰残差同士の相関係数である.

小問[2]は偏相関係数が回帰残差同士の相関係数であることを確かめる問題.

2.2 偏相関係数と偏回帰係数.

単回帰において「単回帰係数=相関係数」であることはよく知られている。
最小二乗法による回帰係数の導出には微積分が必要である。
また、線形代数の行列微分を用いれば重回帰へそのまま拡張できる。
詳細については当サイト回帰分析の流れを参照してほしい。
微積分と線形代数をまだ履修していない場合は結果だけ見ておくとよい。
ここでは、変数が標準化されている場合について確認する。
標準化された回帰分析はパス解析や構造分析につながる重要な手法である。

2.2.1 単回帰.

   \(\begin{pmatrix} y_{1} \\ \vdots \\ y_{n} \end{pmatrix} = \begin{pmatrix} 1 & x_{1}\\ \vdots & \vdots \\ 1 & x_{n} \end{pmatrix} \ \begin{pmatrix} \beta_{1} \\ \beta_{2} \end{pmatrix}+\begin{pmatrix} u_{1} \\ \vdots \\ u_{n} \end{pmatrix}\).

\(X=\begin{pmatrix} 1 & x_{1}\\ \vdots & \vdots \\ 1 & x_{n} \end{pmatrix}\)\(x= \begin{pmatrix} x_{1}\\ \vdots\\ x_{n} \end{pmatrix}\)\(y= \begin{pmatrix} y_{1}\\ \vdots\\ y_{n} \end{pmatrix}\)\(1= \begin{pmatrix} 1\\ \vdots\\ 1 \end{pmatrix}\) とおくと、

   \(y=X\beta +u\)   (行列を用いるとシンプルに表すことができる).

残差 u の大きさが最小となるように回帰係数 \(\beta\) を定めた推定値を \(\hat{\beta}\) とする。

ざっくり表現すれば、
   \(y\fallingdotseq X\hat{\beta}x\) より、 \(X'X\hat{\beta}\fallingdotseq X'y\)   (正規方程式).

<公式>.
   \(\widehat{\beta }=\left( X'X\right) ^{-1}X'y\)  \(\begin{split} Cov\left( \widehat{\beta }\right) =\sigma ^{2}\left( X'X\right) ^{-1} \end{split}\)  (\(\sigma ^{2}\) は誤差分散).

\(X'X=\begin{pmatrix} 1'\\x'\end{pmatrix} \begin{pmatrix} 1 & x \end{pmatrix}= \begin{pmatrix} 1'\cdot 1 & 1'\cdot x \\ x'\cdot1 & x'\cdot x \end{pmatrix} = \begin{pmatrix} n & n\overline x \\ n\overline x & n\overline {x^2} \end{pmatrix}\).

\(X'y=\begin{pmatrix} 1' \\ x'\end{pmatrix}\left( y\right) =\begin{pmatrix} 1'\cdot y \\ x'\cdot y \end{pmatrix} = \begin{pmatrix} n\overline{y} \\ n\overline{xy} \end{pmatrix}\).

   \(\widehat\beta =\begin{pmatrix} \widehat\beta_{1} \\ \widehat\beta_{2} \end{pmatrix}\) とすると、\(\begin{pmatrix} n & n\overline x \\ n\overline x & n\overline {x^2} \end{pmatrix} \begin{pmatrix} \widehat\beta_{1} \\ \widehat\beta_{2} \end{pmatrix} =\begin{pmatrix} n\overline{y} \\ n\overline{xy} \end{pmatrix}\).

よって、
   \(\begin{pmatrix} 1 & \overline{x} \\ \overline{x} & \overline{x^2} \end{pmatrix} \begin{pmatrix} \widehat\beta_{1} \\ \widehat\beta_{2} \end{pmatrix} = \begin{pmatrix} \overline{y} \\ \overline{xy} \end{pmatrix}\)    (正規方程式).

\(\begin{pmatrix} 1 & \overline{x} \\ \overline{x} & \overline{x^2} \end{pmatrix}^{-1} = \dfrac{1}{s_{x}^2} \begin{pmatrix} \overline{x^2} & \ \ -\overline{x} \\ -\overline{x} & \ \ \ \ 1 \end{pmatrix}\) より \(\begin{pmatrix} \widehat\beta_{1} \\ \widehat\beta_{2} \end{pmatrix} =\dfrac{1}{s_{x}^2} \begin{pmatrix} \overline{x^2} & \ \ -\overline{x} \\ -\overline{x} & \ \ \ \ 1 \end{pmatrix} \begin{pmatrix} \overline{y} \\ \overline{xy} \end{pmatrix}\).

\(\begin{pmatrix} \widehat\beta_{1} \\ \widehat\beta_{2} \end{pmatrix} = \dfrac{1}{s_{x}^2} \begin{pmatrix} \overline{x^2}\overline{y}-\overline{x}\cdot \overline{xy} \\ -\overline{x}\overline{y}+\overline{xy} \end{pmatrix} = \dfrac{1}{s_{x}^2} \begin{pmatrix} (s_x^2+\overline{x}^2)\overline{y}-\overline{x}\cdot (s_{xy}+\overline{x}\overline{y})\\ -\overline{x}\overline{y}+\overline{xy} \end{pmatrix} = \dfrac{1}{s_{x}^2} \begin{pmatrix} s_x^2\overline{y}-s_{xy}\cdot \overline{x}\\ s_{xy} \end{pmatrix}\).

  \(\therefore \begin{pmatrix} \widehat\beta_{1} \\ \widehat\beta_{2} \end{pmatrix} = \begin{pmatrix} \overline{y}-\dfrac{s_{xy}}{s_{x}^{2}}\overline{x} \\ \dfrac{s_{xy}}{s_{x}^2} \end{pmatrix} = \begin{pmatrix} \overline{y}-\dfrac{s_{xy}}{s_{x}s_{y}}\dfrac{s_{y}}{s_{x}}\overline{x} \\ \dfrac{s_{xy}}{s_{x}s_{y}}\dfrac{s_{y}}{s_{x}} \end{pmatrix} = \begin{pmatrix} \overline{y}-r_{xy}\dfrac{s_{y}}{s_{x}}\overline{x} \\ r_{xy}\dfrac{s_{y}}{s_{x}} \end{pmatrix}\).

かなり複雑に見えるが、変数が標準化されていれば、もっとシンプルになる。
即ち、
   \(\overline{x^2}\hat{\beta}=\overline{xy}\) から \(\hat{\beta}=\dfrac{\overline{xy}}{x^2}=r_{XY}\).

2.2.2 重回帰.

\(X=\begin{pmatrix} 1 & X_{2}& X_{3} \end{pmatrix}\) は n×3 行列、\(\left( y\right)\) はn×1 行列。
   \(X'X=\begin{pmatrix}1' \\X_{2}'\\X_{3}'\end{pmatrix} \begin{pmatrix} 1 & X_{2}& X_{3} \end{pmatrix} = \begin{pmatrix} 1'\cdot 1 & 1'\cdot X_{2} & 1'\cdot X_{3}\\ X_{2}'\cdot 1 & X_{2}'\cdot X_{2} & X_{2}'\cdot X_{3}\\ X_{3}'\cdot 1& X_{3}'\cdot X_{2} & X_{3}'\cdot X_{3} \end{pmatrix}\).

   \(X'y =\begin{pmatrix} 1' \\ X_{2}'\\ X_{3}' \end{pmatrix}\left( y\right) =\begin{pmatrix} 1'\cdot y \\ X_{2}'\cdot y\\ X_{3}'\cdot y \end{pmatrix}\).

正規方程式は \(\begin{pmatrix} 1'\cdot 1 & 1'\cdot X_{2} & 1'\cdot X_{3}\\ X_{2}'\cdot 1 & X_{2}'\cdot X_{2} & X_{2}'\cdot X_{3}\\ X_{3}'\cdot 1& X_{3}'\cdot X_{2} & X_{3}'\cdot X_{3} \end{pmatrix} \begin{pmatrix} \widehat\beta_{1} \\ \widehat\beta_{2} \\ \widehat\beta_{3} \end{pmatrix} =\begin{pmatrix} 1'\cdot y \\ X_{2}'\cdot y\\ X_{3}'\cdot y \end{pmatrix}\)、即ち.

   \(\begin{pmatrix} n & n\overline{X_{2}} & n\overline{X_{3}}\\ n\overline{X_{2}} &n\overline{X_{2}^2} & n\overline{X_{2}X_{3}}\\ n\overline{X_{3}}& n\overline{X_{2}X_{3}} & n\overline{X_{3}^2} \end{pmatrix} \begin{pmatrix} \widehat\beta_{1} \\ \widehat\beta_{2} \\ \widehat\beta_{3} \end{pmatrix} = \begin{pmatrix} n\overline{y} \\ n\overline{X_{2}y}\\ n\overline{X_{3}y} \end{pmatrix}\)  \(\therefore \begin{pmatrix} 1& \overline{X_{2}} & \overline{X_{3}}\\ \overline{X_{2}} &\overline{X_{2}^2} & \overline{X_{2}X_{3}}\\ \overline{X_{3}}& \overline{X_{2}X_{3}} & \overline{X_{3}^2} \end{pmatrix} \begin{pmatrix} \widehat\beta_{1} \\ \widehat\beta_{2} \\ \widehat\beta_{3} \end{pmatrix} = \begin{pmatrix} \overline{y} \\ \overline{X_{2}y}\\ \overline{X_{3}y} \end{pmatrix}\).

標準化されていれば、正規方程式は相関係数で表される
即ち、
   \(\begin{pmatrix} 1 & r_{23}\\ r_{23}& 1 \end{pmatrix} \begin{pmatrix} \widehat\beta_{2} \\ \widehat\beta_{3} \end{pmatrix} = \begin{pmatrix} r_{2y}\\ r_{3y} \end{pmatrix}\)   \(\therefore \begin{pmatrix} \widehat\beta_{2} \\ \widehat\beta_{3} \end{pmatrix} = \dfrac{1}{1-r_{23}^{2}} \begin{pmatrix} 1 & -r_{23}\\ -r_{23}& 1 \end{pmatrix} \begin{pmatrix} r_{2y}\\ r_{3y} \end{pmatrix}\).

\(\widehat\beta_{2}, \widehat\beta_{3}\) の分散共分散行列は、\(r=cor(X_{2},X_{3})\) として.    \(\dfrac{\sigma ^{2}}{n} \begin{pmatrix} 1 & \ \ \ r\\ r & \ \ \ 1\ \end{pmatrix} ^{-1} = \dfrac{\sigma ^{2}}{n(1-r^2)} \begin{pmatrix} 1 & \ \ \ -r\\ -r & \ \ \ 1\ \end{pmatrix}\).
r=0 ならば、偏回帰係数の分散は \(\dfrac{\sigma^2}{n}\) となる。(データ数が多ければ精度が上がる).

[3]の問題では

   重回帰式は \(Y=b_1X+b_2W\).

変量は標準化されているので、

   \(\begin{pmatrix} b_1 \\ b_2 \end{pmatrix} = \dfrac{1}{1-r_{XW}^{2}} \begin{pmatrix} 1 & -r_{XW}\\ -r_{XW}& 1 \end{pmatrix} \begin{pmatrix} r_{XY}\\ r_{YW} \end{pmatrix}\).

よって、 \(b_1=\dfrac{r_{XY}-r_{XW}r_{YW}}{1-r_{XW}^{2}}\)  \(b_2=\dfrac{r_{YW}-r_{XY}r_{XW}}{1-r_{XW}^{2}}\).

\(r_{XY}=0.72~,~~r_{XW}=0.9~,~~r_{YW}=0.8\) を代入すると、

   \(b_1=\dfrac{r_{XY}-r_{XW}r_{YW}}{1-r_{XW}^{2}}=\dfrac{0.72-0.9\cdot 0.8}{1-0.9^2}=\dfrac{0.152}{0.19}=0.8\).
   \(b_2=\dfrac{r_{YW}-r_{XY}r_{XW}}{1-r_{XW}^{2}}=\dfrac{0.8-0.72\cdot 0.9}{1-0.9^2}=0\).

偏回帰係数と偏相関係数の関係は、
   \(\displaystyle b_1=\dfrac{r_{XY}-r_{XW}r_{YW}}{\sqrt{1-r_{XW}^2}\sqrt{1-r_{YW}^2}} \cdot\sqrt{\frac{1-r_{YW}^2} {1-r_{XW}^2} }= r_{XY/W}\sqrt{\frac{1-r_{YW}^2} {1-r_{XW}^2} }\).
    即ち、偏回帰係数は偏相関係数に共変量の関わりの度合いの調整項を乗じたものになる。

   \(\displaystyle b_1= r_{XY/W}\sqrt{\frac{1-r_{YW}^2} {1-r_{XW}^2} }\)   \(\displaystyle b_2= r_{YW/X}\sqrt{\frac{1-r_{YX}^2} {1-r_{WX}^2} }\).
    まとめ.

  • 相関係数は双方向の回帰係数の相乗平均.
  • 偏相関係数も双方向の標準化偏回帰係数の相乗平均.

[4-1].

\(T_A=X_A+Y_A~,~~cor(X_A,Y_A)=0\) という仮定である。

   \(V(T_A)=V(X_A)+V(Y_A)=2\).

\(Cov(T_A,X_A)=E(T_AX_A)=E(X_A^2+X_AY_A)=1\) であるので、

    \(Cor(T_A,X_A)=\dfrac{Cov(T_A,X_A)}{\sqrt{V(T_A)}\sqrt{V(X_A)}}=\dfrac{1}{\sqrt{2}}\).

同様に、 \(Cor(T_A,Y_A)=\dfrac{Cov(T_A,Y_A)}{\sqrt{V(T_A)}\sqrt{V(Y_A)}}=\dfrac{1}{\sqrt{2}}\).   

[4-2].
解答によると負の相関になるという。
Rでシミュレーションしてもそうなっているが説明できない。
負になる理由は偏相関係数と関係があるかもしれない。

偏相関係数 \(r_{X_AY_A/T_A}=\dfrac{0-1/\sqrt{2}\cdot 1/\sqrt{2}}{\sqrt{1-1/2}\sqrt{1-1/2}}=-\sqrt{1/2}=-0.7\).

自分には説明できないので、Rのシミュレーションにとどめる。

2.3 相関行列からデータ生成.

data<-data.frame(mvrnorm(500, c(0,0,0), matrix(c(1, 0, 1/sqrt(2), 0, 1, 1/sqrt(2), 1/sqrt(2), 1/sqrt(2), 2), 3, 3)))
colnames(data)<-c("x","y","z")
plot(data)

options(rgl.printRglwidget = TRUE) 
library(rgl)
plot3d(data$x, data$y, data$z, type = "s", col = "blue", size = 1)

合計点 z で合否判定をする。 z>2 のデータを抽出する。

subdata <- subset(data, z>1)
plot(subdata$x,subdata$y)

cor(subdata$x,subdata$y)
## [1] -0.2898921

確かに、負の相関となっている。 数理的説明は分からないのでスルーします。

2.4 z=x+y で生成.

データの作り方を変えてみる。

data<-data.frame(mvrnorm(500, c(0,0), matrix(c(1,0,0,1),2,2)))
colnames(data)<-c("x","y")
dat <- mutate(data,z=x+y)
head(dat)
x y z
0.7877469 -0.3443408 0.4434060
-0.3926331 1.2898405 0.8972074
-0.5968431 -0.7223219 -1.3191651
0.1204203 -1.4194864 -1.2990662
0.4183590 0.4727099 0.8910689
2.2013066 -0.7316504 1.4696563
plot3d(dat$x, dat$y, dat$z, type = "s", col = "blue", size = 1)
subdat <- subset(dat, z>1)
plot(subdat$x,subdat$y)

cor(subdat$x,subdat$y)
## [1] -0.6664077

3 問題3_判別分析と2次判別.

sigma <- matrix(c(1,0,0,1), ncol=2)
x1 <- rmvnorm(n=1000, mean=c(1,1), sigma=sigma)  
data1 <- data.frame(X=x1[,1],Y=x1[,2])
x2 <- rmvnorm(n=1000, mean=c(-1,-1), sigma=4*sigma)  
data2 <- data.frame(X=x2[,1],Y=x2[,2])


#model <- lm(Y ~ X, data = data)
plot(x1,asp=1,col="lightblue",xlim=c(-5,5),ylim=c(-5,5))
par(new=T)
plot(x2,asp=1,col="grey88",xlim=c(-5,5),ylim=c(-5,5))

#abline(model)

3.1 (1).

   \(G_1\sim N_2\left(1,{I}\right)\)   \(G_2\sim N_2\left(-1,~2^2{I}\right)\).

   \(D_1^2=(x-1)^T\cdot {I}^{-1}\cdot (x-1)=(x-1)^T (x-1)=|x-1|^2\).

   \(D_2^2=(x+1)^T\cdot ({2^2I})^{-1}\cdot (x+1)=\dfrac{1}{4}(x+1)^T (x+1)=\dfrac{1}{4}|x+1|^2\).

\(x=0\)のとき、
   \(x-1=\begin{pmatrix}-1\\-1\end{pmatrix}\) だから \(D_1^2=2\).
   \(x+1=\begin{pmatrix}1\\1\end{pmatrix}\) だから \(D_2^2=\dfrac{1}{4}\cdot 2=\dfrac{1}{2}\).

3.2 (2).

\(D_1^2=D_2^2\) より、\((x-1)^T (x-1)=\dfrac{1}{4}(x+1)^T (x+1)\).
   \(\begin{split} f(x)&=(x-1)^T (x-1)-\dfrac{1}{4}(x+1)^T (x+1)\\ &=x^Tx-x^T1-1^Tx+1^T1-\dfrac{1}{4}(x^Tx+x^T1+1^Tx+1^T1)\\ &=\dfrac{3}{4}x^Tx-\dfrac{5}{4}(x^T1+1^Tx)+\dfrac{3}{4}1^T1 \end{split}\).

\(x^T1=1^Tx\)(内積)  \(1^T1=1+1=2\) だから.

   \(f(x)=\dfrac{3}{4}x^Tx-\dfrac{5}{2}1^Tx+\dfrac{3}{2}\).

3.3 (3-1).

回転しても分散共分散行列は変わらないから.
回転した座標 \((Y_1,Y_2)\) について、
   \(G_1\sim N_2\left(\begin{pmatrix}\sqrt{2}\\1\end{pmatrix},{I}\right)\)   \(G_2\sim N_2\left(\begin{pmatrix}-\sqrt{2}\\1\end{pmatrix},~2^2{I}\right)\).
G1において.
   \(Y_1\sim N\left(\sqrt{2},~1\right)\).
G2において.
   \(Y_1\sim N\left(-\sqrt{2},~2^2\right)\).

3.4 (3-2).

   \(P(Y_1<0|G_1)=P\left(Z=\dfrac{Y_1-\sqrt{2}}{1}<0\right)=P\left(Z<-\sqrt{2}\right)=0.0793\).
   \(P(Y_1>0|G_2)=P\left(Z=\dfrac{Y_1+\sqrt{2}}{2}>\sqrt{2}\right)=P\left(Z>\dfrac{\sqrt{2}}{2}\right)=0.2389\)     これは、直線 \(Y_1=0\) を判別直線とした線形判別を行った場合の誤判別確率である。

3.4.1 注意.

線形判別分析は、「各群の正規性と等分散性(同じ分散共分散行列を持つ)」を前提としている。

本問では、2群の分散共分散行列が異なっているので2次判別を用いることが適切である。
[3-2]は \(Y_1=0\) を判別直線とした場合の誤判別確率であるが、[4]にあるように、

  • 誤判別確率が大きい.
  • 各群の誤判別確率のバランスが悪い.
    などの問題がある。

3.5 (4).

判別分析は教師あり学習である。学習データを判別すれば、誤判別確率は小さくなると考えられる。従って、テストデータとしては学習データでないものを用いることが望ましい。

これは機械学習の範囲と思う。

4 問題4_回帰分析.

回帰分析.
一見すると単回帰の基本問題に見えるが注意が必要である。
与えられているのは y の予測値 \(\hat{y}\) であるからである。

4.1 回帰の基本.

   予測値  \(\widehat Y=X\widehat\beta\).
   実測値   \(Y=X\widehat \beta+e\)   \(\therefore Y=\widehat Y+e\).
   残差平方和 \(e'e=(Y-X\widehat\beta)'(Y-X\widehat\beta)\).

(全変動)=(回帰変動)+(残差変動)を表す式は、.

   \(\begin{split} (Y-\overline Y\cdot 1)'(Y-\overline Y\cdot 1) &=(\widehat Y-\overline Y\cdot 1)'(\widehat Y-\overline Y\cdot 1)+e'e \end{split}\).

回帰係数を求めるには、\(Cov(x,y)~,~V(y)\) が必要になるが、表で与えられているのは \(\hat{y}\) であることに注意すると、

  • \(Cov(x,y)=Cov(x,\hat{y}+e)=Cov(x,\hat{y})+Cov(x,e)=Cov(x,\hat{y})\).

  • \(V(y)=Cov(\hat{y}+e,\hat{y}+e)=V(\hat{y})+V(e)+2Cov(\hat{y},e)=V(\hat{y})+V(e)\).

4.2 xによる回帰.

Aさんの分析.
<男性>.
   \(E(x)=5(0.1)+6(0.15)+7(0.4)+8(0.35)=7\).
   \(E(x^2)=5^2(0.1)+6^2(0.15)+7^2(0.4)+8^2(0.35)=49.9\).
   \(V(x)=E(x^2)-E(x)^2=0.9\).

   \(E(\hat{y})=7.2(0.1)+7.8(0.15)+8.4(0.4)+9.0(0.35)=8.4\).
   \(E(\hat{y}^2)=7.2^2(0.1)+7.8^2(0.15)+8.4^2(0.4)+9.0^2(0.35)=70.884\).
   \(V(\hat{y})=E(\hat{y}^2)-E(\hat{y})^2=0.324\).
   \(V(y)=S_R+S_e=0.324+0.972=1.296\).

   \(E(x\hat{y})=7.2(0.1)5+7.8(0.15)6+8.4(0.4)7+9.0(0.35)8=59.34\).
   \(Cov(x,\hat{y})=E(x\hat{y})-E(x)E(\hat{y})=59.34-7\times8.4=0.54\).
   \(\therefore ~~Cov(x,y)=Cov(x,\hat{y})=0.54\).

  回帰式を \(y=a+bx\)  とすると、 \(b=\dfrac{Cov(x,y)}{V(x)}~,~~~a=\bar{y}-b\bar{x}\) であるので、
   \(b=\dfrac{Cov(x,y)}{V(x)}=\dfrac{0.54}{0.9}=0.6\).
   \(a=8.4-0.6\times 7=4.2\)   \(\therefore ~~y=4.2+0.6x\).

<女性>.
   \(E(x)=5(0.35)+6(0.4)+7(0.15)+8(0.1)=6\).
   \(E(x^2)=5^2(0.35)+6^2(0.4)+7^2(0.15)+8^2(0.1)=36.9\).
   \(V(x)=E(x^2)-E(x)^2=0.9\).

   \(E(\hat{y})=6.6(0.35)+7.2(0.4)+7.8(0.15)+8.4(0.1)=7.2\).
   \(E(\hat{y}^2)=6.6^2(0.35)+7.2^2(0.4)+7.8^2(0.15)+8.4^2(0.1)=52.164\).
   \(V(\hat{y})=E(\hat{y}^2)-E(\hat{y})^2=0.324\).
   \(V(y)=S_R+S_e=0.324+0.972=1.296\).

   \(E(x\hat{y})=6.6(0.35)5+7.2(0.4)6+7.8(0.15)7+8.4(0.1)8=43.74\).
   \(Cov(x,\hat{y})=E(x\hat{y})-E(x)E(\hat{y})=43.74-6\times7.2=0.54\).
   \(\therefore ~~Cov(x,y)=Cov(x,\hat{y})=0.54\).

  回帰式を \(y=a+bx\)  とすると、 \(b=\dfrac{Cov(x,y)}{V(x)}~,~~~a=\bar{y}-b\bar{x}\) であるので、
   \(b=\dfrac{Cov(x,y)}{V(x)}=\dfrac{0.54}{0.9}=0.6\).
   \(a=7.2-0.6\times 6=3.6\)   \(\therefore ~~y=3.6+0.6x\).

目的は男女に差があるかの検定であり分散分析であるが、
労働時間という共変量があるため、回帰式を比較することになる。

  • 回帰係数の傾きは同じ.
  • 回帰係数の切片が異なる.

これは 共分散分析(ANCOVA)と呼ばれる手法になる。

4.3 1時間あたりの賃金.

Bさんの分析.
  \(\bar{x}_M=7\)   \(\bar{y}_M=8.4\)   \(\bar{y}_M/\bar{x}_M=1.2\).

  \(\bar{x}_F=6\)   \(\bar{y}_F=7.2\)   \(\bar{y}_F/\bar{x}_F=1.2\).

統計量に差は見られない。

4.4 yによる回帰.

<男性>.
 回帰式を \(y=c+dx\)  とすると、 \(d=\dfrac{Cov(x,y)}{V(y)}~,~~~c=\bar{x}-b\bar{y}\) であるので、
   \(d=\dfrac{Cov(x,y)}{V(y)}=\dfrac{0.54}{1.296}=0.4167\).
   \(c=7-0.4167\times 8.4=3.5\)   \(\therefore ~~x=3.5+0.4167y\).

<女性>.
 回帰式を \(y=c+dx\)  とすると、 \(d=\dfrac{Cov(x,y)}{V(y)}~,~~~c=\bar{x}-b\bar{y}\) であるので、
   \(d=\dfrac{Cov(x,y)}{V(y)}=\dfrac{0.54}{1.296}=0.4167\).
   \(c=6-0.4167\times 7.2=3.0\)   \(\therefore ~~x=3.0+0.4167y\).

\(y=7.2\) のとき、

   男性  \(x=3.5+0.4167\times 7.2 =6.5\).
   女性  \(x=3.0+0.4167\times 7.2=6.0\).

同じ賃金を得るのに、男性は 0.5 時間多く働いていることになる。

[4].
解答の通りであるが、あらためて問うまでもないと思う。
出題者が問いたかったのは、

   回帰の残差分散を用いて全分散を算出できるかどうか.

ということのようである。

5 ANCOVAのシミュレーション.

5.1 男性データ.

X=c(5,6,7,8) .

  回帰式:\(Y = 0.6\cdot X+4.2\)   に従ってデータを生成する。

x1 <- rep(c(5,6,7,8),times=c(10,15,40,35))
table(x1)
## x1
##  5  6  7  8 
## 10 15 40 35
e <- rnorm(100,0,sqrt(0.972))
y1 <- 0.6*x1+4.2+e
data <- data.frame(x=x1,y=y1)
plot(data,ylim=c(4,12))

data1 <- cbind(data, z=rep(1,5))
#head(data1)

5.2 女性データ.

X=c(5,6,7,8) .

  回帰式:\(Y = 0.6\cdot X+3.6\)   に従ってデータを生成する。

x2 <- rep(c(5,6,7,8),times=c(35,40,15,10))
table(x2)
## x2
##  5  6  7  8 
## 35 40 15 10
e <- rnorm(100,0,sqrt(0.972))
y2 <- 0.6*x2+3.6+e
data <- data.frame(x=x2,y=y2)
plot(data,ylim=c(4,12))

data2 <- cbind(data, z=rep(2,5))
#head(data1)
data <- rbind(data1,data2)

5.3 傾きの検定.

X1 <- scale(x1,center=T,scale=F)
Y1 <- scale(y1,center=T,scale=F)
X2 <- scale(x2,center=T,scale=F)
Y2 <- scale(y2,center=T,scale=F)

df_sample <- data.frame(group = rep(c("A", "B"), each = 100), y = c(Y1, Y2),
    x = c(X1,X2))
df_sample$group <- df_sample$group %>%
    factor(levels = c("A", "B"))
col=c(1,2)
par(las=1, family="sans", cex=1.2) #図の調整
plot(y ~ x, df_sample, pch=16,
     xlab="x", ylab="y",col=c(1:2))
legend("topleft", legend=c("A", "B"), pch=16, col=1:3, cex=0.8) #凡例の作成
RegressionA <- lm(y ~ x, subset(df_sample, group=="A"))
RegressionB <- lm(y ~ x, subset(df_sample,group=="B"))
Regression <- lm(y ~ x,df_sample)

abline(RegressionA, col=1)
abline(RegressionB, col=2)
abline(Regression, col="lightblue",lwd=5) #全体の回帰

Regression $ coefficients[2]
##         x 
## 0.4553277

5.3.1 平行性検定のF統計量.

s2 <- sum((RegressionA $ residuals)^2)+sum((RegressionB $ residuals)^2)
MSE <- s2/196
MSE
## [1] 1.126785
s1 <- sum(((Regression $ coefficients[2])*(x1-mean(x1))+mean(y1)-y1)^2)+sum(((Regression $ coefficients[2])*(x2-mean(x2))+mean(y2)-y2)^2)
s1/197
## [1] 1.127629
MSA <- (s1-s2)/1
MSA
## [1] 1.292953
f <- MSA/MSE
f
## [1] 1.14747
#F値
1-pf(f, 1, 196)
## [1] 0.2853966

平行性は棄却されない。

library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:psych':
## 
##     logit
## The following object is masked from 'package:dplyr':
## 
##     recode
ANCOVA.m <- lm(y ~ x*group, df_sample) #ANCOVA用線形モデル
Anova(ANCOVA.m)
Sum Sq Df F value Pr(>F)
x 37.318197 1 33.11917 0.0000000
group 0.000000 1 0.00000 1.0000000
x:group 1.292953 1 1.14747 0.2853966
Residuals 220.849932 196 NA NA

5.3.2 切片の検定.

Regression3 <- lm(y ~ x,data)
#summary(Regression3)
s3 <- sum((Regression3 $ residuals)^2)
s3/198
## [1] 1.205562
s3
## [1] 238.7013
s1
## [1] 222.1429
MSA3 <- (s3-s1)/1
MSA3
## [1] 16.55839
f_value <- MSA3/(s1/197)
f_value
## [1] 14.68426
library(car)
ANOVA.m <- lm(y ~ x+factor(z), data) 
Anova(ANOVA.m)
Sum Sq Df F value Pr(>F)
x 37.31820 1 33.09440 0.0000000
factor(z) 16.55839 1 14.68426 0.0001709
Residuals 222.14288 197 NA NA

切片は有意である。