Ch.8 「質的データの統計的分析」の章末問題の解答例 〜自然科学の統計学(東京大学出版会)〜

当記事は「基礎統計学Ⅲ 自然科学の統計学(東京大学出版会)」の読解サポートにあたって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$%であることがわかる。

問題8.5の解答例

まとめ

https://www.amazon.co.jp/dp/4130420674