当記事は「基礎統計学Ⅲ 自然科学の統計学(東京大学出版会)」の読解サポートにあたってChapter.$8$の「質的データの統計的分析」の章末問題の解説について行います。
基本的には書籍の購入者向けの解答例・解説なので、まだ入手されていない方は下記より入手をご検討ください。また、解説はあくまでサイト運営者が独自に作成したものであり、書籍の公式ページではないことにご注意ください。
下記などの演習も並行で抑えておくと良いと思います。
https://www.hello-statisticians.com/practice/stat_practice7.html
章末の演習問題について
問題8.1の解答例
i)
下記を実行することで$\alpha_0, \alpha_1$の推定値の$\hat{\alpha}_0, \hat{\alpha}_1$を得ることができる。
import numpy as np
from scipy import linalg
X = np.ones([5,2])
X[:,1] = np.arange(0.,5.,1)
y = np.array([0.067, 0.267, 0.500, 0.767, 0.900]).reshape([5,1])
X_ = linalg.inv(np.dot(X.T,X))
alpha = np.dot(X_,np.dot(X.T,y))
print(alpha)・実行結果
> print(alpha)
[[ 0.067 ]
[ 0.2166]]上記より、$\hat{\alpha}_0=0.067, \hat{\alpha}_1=0.2166$が得られる。
ⅱ)
下記を実行することで結果を得ることができる。
import numpy as np
from scipy import linalg
X = np.ones([5,2])
X[:,1] = np.arange(0.,5.,1)
y = np.array([0.067, 0.267, 0.500, 0.767, 0.900]).reshape([5,1])
X_ = linalg.inv(np.dot(X.T,X))
alpha = np.dot(X_,np.dot(X.T,y))
print(np.dot(X,alpha))・実行結果
> print(np.dot(X,alpha))
[[ 0.067 ]
[ 0.2836]
[ 0.5002]
[ 0.7168]
[ 0.9334]]ⅲ)
省略
iv)
$0.2166X+0.067>1, 0.2166X+0.067<0$をそれぞれ解けばよい。
・$0.2166X+0.067>1$
$X>4.3074…$
・$0.2166X+0.067<0$
$X<-0.3093…$
問題8.2の解答例
i)
$\displaystyle P = F(\beta_{0}+\beta_{1}X) = \frac{1}{2}$は$\beta_{0}+\beta_{1}X=0$の時に成立する。よってこのときの$X$の値は下記のように導出できる。
$$
\large
\begin{align}
\beta_{0} + \beta_{1}X &= 0 \\
\beta_{1}X &= -\beta_{0} \\
X &= – \frac{\beta_{0}}{\beta_{1}}
\end{align}
$$
ⅱ)
それぞれ下記のように考えることができる。
・プロビット
$$
\large
\begin{align}
P &= F(\beta_{0}+\beta_{1}X) \\
F(Z) &= \int_{-\infty}^{Z} \frac{1}{\sqrt{2 \pi}} e^{-\frac{x^2}{2}} dx
\end{align}
$$
上記を元に$\displaystyle \frac{dP}{dX}$を計算する。
$$
\large
\begin{align}
\frac{dP}{dX} &= \frac{dP}{dZ} \cdot \frac{dZ}{dX} \\
&= \frac{1}{\sqrt{2 \pi}} e^{-\frac{Z^2}{2}} \cdot \beta_{1}
\end{align}
$$
上記を$Z$の関数と見ると、$Z=0$のとき最大値を取ることがわかる。よって、$\displaystyle \left| \frac{dP}{dX} \right|$は$\displaystyle Z=0, X = – \frac{\beta_{0}}{\beta_{1}}$のとき下記のような最大値を取る。
$$
\large
\begin{align}
\frac{1}{\sqrt{2 \pi}} e^{-\frac{0^2}{2}} \cdot \beta_{1} &= \frac{1}{\sqrt{2 \pi}} \cdot \beta_{1} \\
&= 0.3989 \beta_{1}
\end{align}
$$
・ロジット
$$
\large
\begin{align}
P &= F(\beta_{0}+\beta_{1}X) \\
F(Z) &= \frac{e^{Z}}{1+e^{Z}} \\
&= \frac{1}{1+e^{-Z}}
\end{align}
$$
$Z=\beta_{0}+\beta_{1}X$とおき、上記を元に$\displaystyle \frac{dP}{dX}$を計算する。
$$
\large
\begin{align}
\frac{dP}{dX} &= \frac{dP}{dZ} \cdot \frac{dZ}{dX} \\
&= -\frac{1}{(1+e^{-Z})^2} \cdot -e^{-Z} \cdot \beta_{1} \\
&= \frac{e^{-Z}}{(1+e^{-Z})^2} \cdot \beta_{1}
\end{align}
$$
ここで上記を最大にする$Z$がわかれば$Z=\beta_{0}+\beta_{1}X$より対応する$X$の値も計算できる。よって、以下では$\displaystyle g(Z) = \frac{e^{-Z}}{(1+e^{-Z})^2}$を最大にする$Z$を考える。
$$
\large
\begin{align}
\frac{dg(Z)}{dZ} &= \frac{d}{dZ} ((1+e^{-Z})^{-2}e^{-Z}) \\
&= -2(1+e^{-Z})^{-3} \times -e^{-Z} \times e^{-Z} + (1+e^{-Z})^{-2} \times -e^{-Z} \\
&= \frac{e^{-Z}(e^{-Z}-1)}{2(1+e^{-Z})^{3}}
\end{align}
$$
上記は単調減少かつ$Z=0$で$\displaystyle \frac{dg(Z)}{dZ}=0$なので$Z=0$のとき$g(Z)$は最大値を取る。このことは$\displaystyle \left| \frac{dP}{dX} \right|$は$\displaystyle Z=0, X = – \frac{\beta_{0}}{\beta_{1}}$のとき下記のような最大値を取ることを意味する。
$$
\large
\begin{align}
\frac{e^{-0}}{(1+e^{-0})^2} \cdot \beta_{1} &= \frac{1}{(1+1)^2} \cdot \beta_{1} \\
&= 0.25 \beta_{1}
\end{align}
$$
問題8.3の解答例
例8.2と例8.5に関してそれぞれ確認を行う。
・例8.2
print(2.324/1.086)
print(1.162/0.543)・実行結果
> print(2.324/1.086)
2.13996316759
> print(1.162/0.543)
2.13996316759・例8.5
print(50.26/27.56)
print(2.526/1.358)
print(0.612/0.337)・実行結果
> print(50.26/27.56)
1.8236574746
> print(2.526/1.358)
1.86008836524
> print(0.612/0.337)
1.81602373887$\displaystyle \frac{\pi}{\sqrt{3}} \simeq 1.81$より、例8.5で概ね成立していることが確認できる。
問題8.4の解答例
下記を実行することで結果を得ることができる。
import numpy as np
X = np.array([[0.,89.,0.], [0.,85.,0.], [1.,84.,1.], [0.,83.,0.], [1.,81.,1.], [0.,88.,1.], [1.,87.,0.], [1.,86.,0.], [0.,86.,1.], [1.,84.,10.], [0.,87.,0.], [1.,83.,1.], [1.,82.,1.], [0.,83.,0.], [0.,82.,0.], [0.,86.,0.], [0.,78.,0.], [0.,77.,1.], [0.,85.,1.], [1.,75.,1.], [1.,80.,0.], [1.,82.,0.], [1.,79.,1.], [0.,83.,1.], [0.,82.,1.], [0.,89.,0.], [1.,86.,1.], [1.,80.,1.], [1.,81.,0.], [0.,78.,0.]])
y = np.array([1., 1., 0., 1., 0., 1., 0., 0., 1., 0., 1., 0., 0., 1., 0., 1., 0., 0., 1., 0., 0., 0., 0., 0., 0., 1., 1., 0., 1., 0.])
y_probit = -27.56 - 1.358*X[:,0] + 0.337*X[:,1] - 0.447*X[:,2]
y_probit[y_probit>0] = 1
y_probit[y_probit<=0] = 0
y_logi = -50.26 - 2.526*X[:,0] + 0.612*X[:,1] - 0.610*X[:,2]
y_logi[y_logi>0] = 1
y_logi[y_logi<=0] = 0
print(np.float(np.sum(y==y_probit))/y.shape[0])
print(np.float(np.sum(y==y_logi))/y.shape[0])・実行結果
> print(np.float(np.sum(y==y_probit))/y.shape[0])
0.833333333333
> print(np.float(np.sum(y==y_logi))/y.shape[0])
0.9上記より正答率はプロビットに対して$83.3$％、ロジットに対して$90.0$％であることがわかる。