ブログ

「統計学実践ワークブック」 演習問題etc Ch.26 「その他の多変量解析手法」

当記事は「統計学実践ワークブック(学術図書出版社)」の読解サポートにあたってChapter.26の「その他の多変量解析手法」に関して演習問題を中心に解説を行います。多次元尺度構成法や数量化法、対応分析は用いられることがあるので、演習を通して抑えておくと良いと思われました。

本章のまとめ

多次元尺度構成法

下記で詳しく取り扱った。
https://www.hello-statisticians.com/explain-terms-cat/mds1.html

演習問題解説

問26.1

$[1]$
$$
\large
\begin{align}
B = – \frac{1}{2} \left( \mathit{I}_{n} – \frac{1}{n} \mathit{J}_n \right) D \left( \mathit{I}_{n} – \frac{1}{n} \mathit{J}_n \right)
\end{align}
$$

上記の式に基づいて$B$の計算を行う。

import numpy as np

I_4 = np.eye(4)
J_4 = np.ones([4,4])
D = np.array([[0., 5., 5., 16.], [5., 0., 4., 5.], [5., 4., 0., 5.], [16., 5., 5., 0.]])

B = -np.dot(np.dot(I_4-J_4/4., D), I_4-J_4/4.)/2.
lamb, eig_vec = linalg.eig(B)

print("B:")
print(B)
print("lambda: {}".format(lamb.real))
print("eigen vector:")
print(eig_vec.real)

・実行結果

B:
[[ 4. -0. -0. -4.]
 [-0.  1. -1. -0.]
 [-0. -1.  1. -0.]
 [-4. -0. -0.  4.]]
lambda: [ 8.  0.  2.  0.]
eigen vector:
array([[ 0.70710678, -0.70710678,  0.        ,  0.        ],
       [ 0.        ,  0.        , -0.70710678, -0.70710678],
       [ 0.        ,  0.        ,  0.70710678, -0.70710678],
       [-0.70710678, -0.70710678,  0.        ,  0.        ]])

$[2]$

下記のようにそれぞれ計算できる。

import numpy as np
from scipy import linalg

I_4 = np.eye(4)
J_4 = np.ones([4,4])
D = np.array([[0., 5., 5., 16.], [5., 0., 4., 5.], [5., 4., 0., 5.], [16., 5., 5., 0.]])

B = -np.dot(np.dot(I_4-J_4/4., D), I_4-J_4/4.)/2.
lamb, eig_vec = linalg.eig(B)

x = eig_vec.real[:,0]*np.sqrt(lamb.real[0])
x_row = np.repeat(x,4).reshape([4,4])
x_col = np.repeat(x,4).reshape([4,4]).T
D = (x_row-x_col)**2

print("each x: {}".format(x))
print("D:")
print(D)

・実行結果

each x: [ 2.  0.  0. -2.]
D:
[[  0.   4.   4.  16.]
 [  4.   0.   0.   4.]
 [  4.   0.   0.   4.]
 [ 16.   4.   4.   0.]]

上記で計算したDでは元々の$D$を完全には再現されていないことが確認できる。

$[3]$
下記のようにそれぞれ計算できる。

import numpy as np
from scipy import linalg

I_4 = np.eye(4)
J_4 = np.ones([4,4])
D = np.array([[0., 5., 5., 16.], [5., 0., 4., 5.], [5., 4., 0., 5.], [16., 5., 5., 0.]])

B = -np.dot(np.dot(I_4-J_4/4., D), I_4-J_4/4.)/2.
lamb, eig_vec = linalg.eig(B)

x = np.dot(eig_vec.real[:,np.array([0,2])], np.array([[np.sqrt(lamb.real[0]), 0.], [0., np.sqrt(lamb.real[2])]]))
D = np.zeros([4,4])
for i in range(4):
    for j in range(4):
        D[i,j] = np.sum((x[i,:]-x[j,:])**2)

print("each x:")
print(x)
print("D:")
print(D)

・実行結果

each x:
[[ 2.  0.]
 [ 0. -1.]
 [ 0.  1.]
 [-2.  0.]]
D:
[[  0.   5.   5.  16.]
 [  5.   0.   4.   5.]
 [  5.   4.   0.   5.]
 [ 16.   5.   5.   0.]]

xの$2$行目と$3$行目がワークブックの解答とは逆だが、固有ベクトルの向きが逆であるだけであり、該当する指定が特にないので、上記の結果も正しいと考えることができる。また、上記で計算したDでは元々の$D$が完全に再現されていることも確認できる。

参考

・準1級関連まとめ
https://www.hello-statisticians.com/toukeikentei-semi1

「統計学実践ワークブック」 演習問題etc Ch.23 「判別分析」

当記事は「統計学実践ワークブック(学術図書出版社)」の読解サポートにあたってChapter.$23$の「判別分析」に関して演習問題を中心に解説を行います。フィッシャーの判別分析やSVMなどはよく出てくる手法なので、演習を通して抑えておくと良いと思われました。

本章のまとめ

フィッシャーの線形判別の導出

フィッシャーの線形判別」で詳しく取り扱った。

$2$次判別分析

「フィッシャーの線形判別」などの線形判別では$2$クラスの共分散行列が等しいと仮定するが、この際にそれぞれのクラスの共分散行列を考えることで$2$次判別分析となる。

以下、サンプルを$\mathbf{x} \in \mathbb{R}^{2}$、クラス$1$を$y=1$、クラス$2$を$y=-1$、クラス$1$の平均を$\mu_{1} \in \mathbb{R}^{2}$共分散行列を$\Sigma_1 \in \mathbb{R}^{2 \times 2}$、クラス$2$の平均を$\mu_{2} \in \mathbb{R}^{2}$共分散行列を$\Sigma_2 \in \mathbb{R}^{2 \times 2}$でそれぞれ表すことを考える。

このとき、判別関数$f(\mathbf{x})$を下記のように定義する。
$$
\large
\begin{align}
f(\mathbf{x}) &= \log{ \frac{P(y=1|\mathbf{x})}{P(y=-1|\mathbf{x})} } \\
&= \log{ \frac{P(\mathbf{x}|y=1)P(y=1)/\cancel{P(\mathbf{x})}}{P(\mathbf{x}|y=-1)P(y=-1)/\cancel{P(\mathbf{x})}} } \\
&= \log{ \frac{P(\mathbf{x}|y=1)P(y=1)}{P(\mathbf{x}|y=-1)P(y=-1)} } \quad (1)
\end{align}
$$

ここで上記の$(1)$式に対し、下記が成立すると仮定する。
$$
\large
\begin{align}
P(\mathbf{x}|y=1) &= \mathcal{N}(\mathbf{x}|\mu_{1},\Sigma_{1}) \\
&= \frac{1}{(2 \pi)^{2/2} |\Sigma_{1}|^{1/2}} \exp \left( (\mathbf{x}-\mu_{1})^{\mathrm{T}}\Sigma_{1}^{-1}(\mathbf{x}-\mu_{1}) \right) \\
&= \frac{1}{2 \pi |\Sigma_{1}|^{1/2}} \exp \left( (\mathbf{x}-\mu_{1})^{\mathrm{T}}\Sigma_{1}^{-1}(\mathbf{x}-\mu_{1}) \right) \\
P(\mathbf{x}|y=-1) &= \mathcal{N}(\mathbf{x}|\mu_{2},\Sigma_{2}) \\
&= \frac{1}{(2 \pi)^{2/2} |\Sigma_{2}|^{1/2}} \exp \left( (\mathbf{x}-\mu_{1})^{\mathrm{T}}\Sigma_{2}^{-1}(\mathbf{x}-\mu_{2}) \right) \\
&= \frac{1}{2 \pi |\Sigma_{2}|^{1/2}} \exp \left( (\mathbf{x}-\mu_{1})^{\mathrm{T}}\Sigma_{2}^{-1}(\mathbf{x}-\mu_{2}) \right) \\
P(y=1) &= \pi_{1} \\
P(y=-1) &= \pi_{2}
\end{align}
$$

このとき、$(1)$式は下記のように変形を行える。
$$
\large
\begin{align}
& f(\mathbf{x}) = \log{\frac{P(y=1|\mathbf{x})}{P(y=-1|\mathbf{x})}} \quad (1) \\
&= \log{\frac{P(\mathbf{x}|y=1)P(y=1)/\cancel{P(\mathbf{x})}}{P(\mathbf{x}|y=-1)P(y=-1)/\cancel{P(\mathbf{x})}}} \\
&= \log{\frac{\mathcal{N}(\mathbf{x}|\mu_1,\Sigma_1)\pi_{1}}{\mathcal{N}(\mathbf{x}|\mu_2,\Sigma_2)\pi_{2}}} \\
&= \log{\frac{\pi_{1}}{\pi_{2}}} + \log{ \left[ \frac{\cancel{2 \pi}}{\cancel{2 \pi}}\frac{|\Sigma_2|^{1/2}}{|\Sigma_1|^{1/2}}\exp \left( -\frac{1}{2}(\mathbf{x}-\mu_1)^{\mathrm{T}}\Sigma_1^{-1}(\mathbf{x}-\mu_1) +\frac{1}{2}(\mathbf{x}-\mu_2)^{\mathrm{T}}\Sigma_2^{-1}(\mathbf{x}-\mu_2) \right) \right] } \\
&= -\frac{1}{2}(\mathbf{x}-\mu_1)^{\mathrm{T}}\Sigma_1^{-1}(\mathbf{x}-\mu_1) + \frac{1}{2}(\mathbf{x}-\mu_2)^{\mathrm{T}}\Sigma_2^{-1}(\mathbf{x}-\mu_2) + \log{\frac{\pi_{1}}{\pi_{2}}} + \log{\frac{|\Sigma_2|^{1/2}}{|\Sigma_1|^{1/2}}} \\
&= \mathbf{x}^{\mathrm{T}} \frac{1}{2}(-\Sigma_1^{-1}+\Sigma_2^{-1}) \mathbf{x} + (\mu_1^{\mathrm{T}}\Sigma_1^{-1}-\mu_2^{\mathrm{T}}\Sigma_2^{-1}) \mathbf{x} + \frac{1}{2}(-\mu_1^{\mathrm{T}}\Sigma_1^{-1}\mu_1+\mu_2^{\mathrm{T}}\Sigma_2^{-1}\mu_2) \\
& \qquad + \log{\frac{\pi_{1}}{\pi_{2}}} + \log{\frac{|\Sigma_2|^{1/2}}{|\Sigma_1|^{1/2}}} \quad (2)
\end{align}
$$

上記$(2)$式の$2$次の項が$\displaystyle \mathbf{x}^{\mathrm{T}} \frac{1}{2}(-\Sigma_1^{-1}+\Sigma_2^{-1}) \mathbf{x}$で表されることから、$\Sigma_1 \neq \Sigma_2$のとき、判別関数$f(\mathbf{x})$が$2$次であることが確認できる。

一方で$\Sigma_1 = \Sigma_2$が成立するとき、$(2)$式は$1$次式となるので、ここでの判別は線形判別に対応すると考えられる。

サポートベクトルマシン

演習問題解説

問23.1

$[1]$
$2$つの群が直線で分割できない際はフィッシャーの線形判別よりもガウシアンカーネルを用いたSVMの方が良い判別精度を示すと考えることができる。

$[2]$
線形判別の方がテストケースに対して良い結果を示していることから、過学習が生じていることが読み取れる。この対策にあたってはクロスバリデーションなどを行うことで、パラメータを適切な値に調節することなどが考えられる。

問23.2

下記を実行することで分類を行うことができる。

import numpy as np

w1 = np.array([-0.1, -0.1, 0.2, 0.4, 0.9])
w2 = np.array([0.0, 0.2, -0.2, 0.9, -0.3])

ave_x1 = np.array([1.2, 1.4, 0.3, 6.0, 4.6])
ave_x2 = np.array([1.4, 1.4, 0.3, 0.7, 2.1])
ave_x3 = np.array([0.8, 0.3, 1.6, 3.6, 9.2])

vec_x1 = np.array([np.dot(w1,ave_x1), np.dot(w2,ave_x1)])
vec_x2 = np.array([np.dot(w1,ave_x2), np.dot(w2,ave_x2)])
vec_x3 = np.array([np.dot(w1,ave_x3), np.dot(w2,ave_x3)])

new_x = np.array([2.5, -1.1, 0.0, -0.8, 1.0])
vec_new_x = np.array([np.dot(w1,new_x), np.dot(w2,new_x)])

cat_class = ["E", "F", "Q"]

dist2_vec = np.array([np.sum((vec_new_x-vec_x1)**2), np.sum((vec_new_x-vec_x2)**2), np.sum((vec_new_x-vec_x3)**2)])

print("Class: {}".format(cat_class[np.argmin(dist2_vec)]))

・実行結果

> print("Class: {}".format(cat_class[np.argmin(dist2_vec)]))
Class: F

問23.3

$[1]$
$$
\large
\begin{align}
f_{q}(\mathbf{x}) = \log{\frac{P(y=1|\mathbf{x})}{P(y=-1|\mathbf{x})}} \quad (1)
\end{align}
$$

ここで$P(\mathbf{x}|y=1)=\mathcal{N}(\mathbf{x}|\mu_1,\Sigma_1), P(\mathbf{x}|y=-1)=\mathcal{N}(\mathbf{x}|\mu_2,\Sigma_2)$より、$(1)$式は下記のように変形できる。
$$
\large
\begin{align}
& f_{q}(\mathbf{x}) = \log{\frac{P(y=1|\mathbf{x})}{P(y=-1|\mathbf{x})}} \quad (1) \\
&= \log{\frac{P(\mathbf{x}|y=1)P(y=1)/\cancel{P(\mathbf{x})}}{P(\mathbf{x}|y=-1)P(y=-1)/\cancel{P(\mathbf{x})}}} \\
&= \log{\frac{\mathcal{N}(\mathbf{x}|\mu_1,\Sigma_1)\pi_{1}}{\mathcal{N}(\mathbf{x}|\mu_2,\Sigma_2)\pi_{2}}} \\
&= \log{\frac{\pi_{1}}{\pi_{2}}} + \log{ \left[ \frac{\cancel{2 \pi}}{\cancel{2 \pi}}\frac{|\Sigma_2|^{1/2}}{|\Sigma_1|^{1/2}}\exp \left( -\frac{1}{2}(\mathbf{x}-\mu_1)^{\mathrm{T}}\Sigma_1^{-1}(\mathbf{x}-\mu_1) + \frac{1}{2}(\mathbf{x}-\mu_2)^{\mathrm{T}}\Sigma_2^{-1}(\mathbf{x}-\mu_2) \right) \right] } \\
&= -\frac{1}{2}(\mathbf{x}-\mu_1)^{\mathrm{T}}\Sigma_1^{-1}(\mathbf{x}-\mu_1) + \frac{1}{2}(\mathbf{x}-\mu_2)^{\mathrm{T}}\Sigma_2^{-1}(\mathbf{x}-\mu_2) + \log{\frac{\pi_{1}}{\pi_{2}}} + \log{\frac{|\Sigma_2|^{1/2}}{|\Sigma_1|^{1/2}}} \\
&= \mathbf{x}^{\mathrm{T}} \frac{1}{2}(-\Sigma_1^{-1}+\Sigma_2^{-1}) \mathbf{x} + (\mu_1^{\mathrm{T}}\Sigma_1^{-1}-\mu_2^{\mathrm{T}}\Sigma_2^{-1}) \mathbf{x} + \frac{1}{2}(-\mu_1^{\mathrm{T}}\Sigma_1^{-1}\mu_1+\mu_2^{\mathrm{T}}\Sigma_2^{-1}\mu_2) \\
& \qquad + \log{\frac{\pi_{1}}{\pi_{2}}} + \log{\frac{|\Sigma_2|^{1/2}}{|\Sigma_1|^{1/2}}}
\end{align}
$$

上記より$A, \mathbf{b}, c$はそれぞれ下記のように得られる。
$$
\large
\begin{align}
A &= \frac{1}{2}(-\Sigma_1^{-1}+\Sigma_2^{-1}) \\
\mathbf{b} &= \Sigma_1^{-1}\mu_1 – \Sigma_2^{-1}\mu_2 \\
c &= \frac{1}{2}(-\mu_1^{\mathrm{T}}\Sigma_1^{-1}\mu_1+\mu_2^{\mathrm{T}}\Sigma_2^{-1}\mu_2) + \log{\frac{\pi_{1}}{\pi_{2}}} + \log{\frac{|\Sigma_2|^{1/2}}{|\Sigma_1|^{1/2}}}
\end{align}
$$

$[2]$
与えられた$\Sigma_1,\Sigma_2$を元に$\Sigma_1^{-1},\Sigma_2^{-1}$はそれぞれ下記のように得られる。
$$
\large
\begin{align}
\Sigma_1^{-1} &= \left(\begin{array}{cc} 0.8 & 0 \\ 0 & 0.5 \end{array} \right)^{-1} \\
&= \left(\begin{array}{cc} 1/0.8 & 0 \\ 0 & 1/0.5 \end{array} \right) \\
&= \left(\begin{array}{cc} 1.25 & 0 \\ 0 & 2 \end{array} \right) \\
\Sigma_2^{-1} &= \left(\begin{array}{cc} 0.4 & 0 \\ 0 & 0.25 \end{array} \right)^{-1} \\
&= \left(\begin{array}{cc} 1/0.4 & 0 \\ 0 & 1/0.25 \end{array} \right) \\
&= \left(\begin{array}{cc} 2.5 & 0 \\ 0 & 4 \end{array} \right)
\end{align}
$$

また、$\pi_{1}=100/300=1/3, \pi_{2}=200/300=2/3$が成立する。その他与えられた数値を用いることで、$[1]$の$A, \mathbf{b}, c$はそれぞれ下記のように計算できる。

$$
\large
\begin{align}
A &= \frac{1}{2}(-\Sigma_1^{-1}+\Sigma_2^{-1}) \\
&= \frac{1}{2} \left( -\left(\begin{array}{cc} 1.25 & 0 \\ 0 & 2 \end{array} \right)+\left(\begin{array}{cc} 2.5 & 0 \\ 0 & 4 \end{array} \right) \right) \\
&= \frac{1}{2} \left(\begin{array}{cc} 1.25 & 0 \\ 0 & 2 \end{array} \right) = \left(\begin{array}{cc} 0.625 & 0 \\ 0 & 1 \end{array} \right) \\
\mathbf{b} &= \Sigma_1^{-1}\mu_1 – \Sigma_2^{-1}\mu_2 \\
&= \left(\begin{array}{cc} 1.25 & 0 \\ 0 & 2 \end{array} \right)\left(\begin{array}{c} 2 \\ 2 \end{array} \right) – \left(\begin{array}{cc} 2.5 & 0 \\ 0 & 4 \end{array} \right)\left(\begin{array}{c} 0 \\ 0 \end{array} \right) \\
&= \left(\begin{array}{c} 2.5 \\ 4 \end{array} \right) \\
c &= \frac{1}{2}(-\mu_1^{\mathrm{T}}\Sigma_1^{-1}\mu_1+\mu_2^{\mathrm{T}}\Sigma_2^{-1}\mu_2) + \log{\frac{\pi_{1}}{\pi_{2}}} + \log{\frac{|\Sigma_2|^{1/2}}{|\Sigma_1|^{1/2}}} \\
&= – \frac{1}{2}\left(\begin{array}{cc} 2 & 2 \end{array} \right)\left(\begin{array}{cc} 1.25 & 0 \\ 0 & 2 \end{array} \right)\left(\begin{array}{c} 2 \\ 2 \end{array} \right) + \log{\frac{1/3}{2/3}} + \log{\frac{0.1^{1/2}}{0.4^{1/2}}} \\
&= – \left(\begin{array}{cc} 1 & 1 \end{array} \right)\left(\begin{array}{c} 2.5 \\ 4 \end{array} \right) + \log{\frac{1}{2}} + \frac{1}{2}\log{\frac{1}{4}} \\
&= – 2.5 – 4 – \log{2} – \log{2} \\
&= -7.886…
\end{align}
$$

よって、判別関数を$f(\mathbf{x})$とおくと、$f(\mathbf{x})$は下記のように表せる。
$$
\large
\begin{align}
f(\mathbf{x}) &= \mathbf{x}^{\mathrm{T}} \left(\begin{array}{cc} 0.625 & 0 \\ 0 & 1 \end{array} \right) \mathbf{x} + \left(\begin{array}{c} 2.5 \\ 4 \end{array} \right)^{\mathrm{T}} \mathbf{x} – 7.886… \\
&= \mathbf{x}^{\mathrm{T}} \left(\begin{array}{cc} 0.625 & 0 \\ 0 & 1 \end{array} \right) \mathbf{x} + \left(\begin{array}{cc} 2.5 & 4 \end{array} \right) \mathbf{x} – 7.886…
\end{align}
$$

上記に対して$\mathbf{x} = \left(\begin{array}{cc} 1 & 1 \end{array} \right)^{\mathrm{T}}, \left(\begin{array}{cc} 0 & 1 \end{array} \right)^{\mathrm{T}}$を代入すると下記のように計算できる。
・$\mathbf{x} = \left(\begin{array}{cc} 1 & 1 \end{array} \right)^{\mathrm{T}}$
$$
\large
\begin{align}
f(\mathbf{x}) &= \mathbf{x}^{\mathrm{T}} \left(\begin{array}{cc} 0.625 & 0 \\ 0 & 1 \end{array} \right) \mathbf{x} + \left(\begin{array}{cc} 2.5 & 4 \end{array} \right) \mathbf{x} – 7.886… \\
&= \left(\begin{array}{cc} 1 & 1 \end{array} \right) \left(\begin{array}{cc} 0.625 & 0 \\ 0 & 1 \end{array} \right) \left(\begin{array}{c} 1 \\ 1 \end{array} \right) + \left(\begin{array}{cc} 2.5 & 4 \end{array} \right) \left(\begin{array}{c} 1 \\ 1 \end{array} \right) – 7.886… \\
&= 1+0.625+2.5+4-7.886… = 0.2387… > 0
\end{align}
$$

よってクラス$1$に分類される。

・$\mathbf{x} = \left(\begin{array}{cc} 0 & 1 \end{array} \right)^{\mathrm{T}}$
$$
\large
\begin{align}
f(\mathbf{x}) &= \mathbf{x}^{\mathrm{T}} \left(\begin{array}{cc} 0.625 & 0 \\ 0 & 1 \end{array} \right) \mathbf{x} + \left(\begin{array}{cc} 2.5 & 4 \end{array} \right) \mathbf{x} – 7.886… \\
&= \left(\begin{array}{cc} 0 & 1 \end{array} \right) \left(\begin{array}{cc} 0.625 & 0 \\ 0 & 1 \end{array} \right) \left(\begin{array}{c} 0 \\ 1 \end{array} \right) + \left(\begin{array}{cc} 2.5 & 4 \end{array} \right) \left(\begin{array}{c} 0 \\ 1 \end{array} \right) – 7.886… \\
&= 1+4-7.886… = -2.886 < 0
\end{align}
$$

よってクラス$2$に分類される。

参考

・準$1$級関連まとめ
https://www.hello-statisticians.com/toukeikentei-semi1

「統計学実践ワークブック」 演習問題etc Ch.24 「クラスター分析」

当記事は「統計学実践ワークブック(学術図書出版社)」の読解サポートにあたってChapter.24の「クラスター分析」に関して演習問題を中心に解説を行います。実際に分析などを行うにあたってよく出てくる手法なので、演習を通して抑えておくと良いと思われました。

本章のまとめ

演習問題解説

問24.1

$[1]$
表$24.2$を元に考えると、$3,4$が最も近く、$3,4$と$5$が最近隣法でも最遠隣法でも最も近い。よって、$3,4$に$A, C$が対応し、$5$に$F$が対応することがわかる。次に相互に近いのが$1-2$であり、これが$D, E$に対応することがわかる。また、$3,4,5$と$6$間は最近隣法だと$2.7$、最遠隣法だと$4.6$が対応する一方で、$3,4,5$と$1,2$間は最近隣法だと$2.5$、最遠隣法だと$5.6$が対応する。よって、最近隣法は$(b)$、最遠隣法は$(a)$が対応する。

また、番号とアルファベットの対応は下記であることもわかる。

1,2: D,E
3,4: A,C
5:  F
6:  B

$[2]$
下記を実行することで計算することができる。

import numpy as np

X = np.array([[7.7, 5.8, 7.7], [10.1, 5.9, 7.3], [4.9, 4.1, 7.9], [4.8, 5.0, 9.0], [6.4, 3.7, 8.0], [4.0, 5.2, 11.6]])
x1_average0 = np.array([9., 6., 7.])
x2_average0 = np.array([5., 5., 10.])

sample_x1 = np.sum((X-x1_average0)**2,axis=1)>=np.sum((X-x2_average0)**2,axis=1)
sample_x2 = np.sum((X-x1_average0)**2,axis=1)<np.sum((X-x2_average0)**2,axis=1)

x1_average1 = np.mean(X[sample_x1,:],axis=0)
x2_average1 = np.mean(X[sample_x2,:],axis=0)

print("average x1: {}".format(x1_average1))
print("average x2: {}".format(x2_average1))

・実行結果

> print("average x1: {}".format(x1_average1))
average x1: [ 5.025  4.5    9.125]
> print("average x2: {}".format(x2_average1))
average x2: [ 8.9   5.85  7.5 ]

参考

・準1級関連まとめ
https://www.hello-statisticians.com/toukeikentei-semi1

Pythonを用いた確率過程のプログラミング 〜ブラウン運動、ポアソン過程 etc〜

数式だけの解説ではわかりにくい場合もあると思われるので、統計学の手法や関連する概念をPythonのプログラミングで表現します。当記事では確率過程(stochastic process)の理解にあたって、ブラウン運動やポアソン過程のPythonでの実装を取り扱いました。

・Pythonを用いた統計学のプログラミングまとめ
https://www.hello-statisticians.com/stat_program

ブラウン運動

標準ブラウン運動・ウィーナー過程

標準ブラウン運動$N(0, t)$はウィーナー過程と呼ばれることもある。標準ブラウン運動は下記のように表現できる。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(0)

x = np.arange(0,100,1)
y = np.zeros([x.shape[0],5])

y[0,:] = stats.norm.rvs(size=5)
for i in range(x.shape[0]-1):
    y[i+1,:] = y[i,:] + stats.norm.rvs(size=5)

for i in range(y.shape[1]):
    plt.plot(x,y[:,i])

plt.show()

・実行結果

正規乱数の生成を簡易化するにあたってscipy.stats.norm.rvsを用いた。パッケージを用いない際は乗算合同法やM系列、メルセンヌツイスタ法などを用いて作成した一様乱数を元に中心極限定理の応用やボックス・ミュラー法を用いて生成することができる。
https://www.hello-statisticians.com/explain-terms-cat/random_sampling1.html
https://www.hello-statisticians.com/explain-books-cat/toukeigaku-aohon/stat_blue_ch11.html

ブラウン運動

確率過程の生成

下記を実行することでブラウン運動$N(\mu t, t)$を生成することができる。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(0)

mu, sigma = 0.2, 1
x = np.arange(0,100,1)
y = np.zeros([x.shape[0],5])

y[0,:] = stats.norm.rvs(size=5)
for i in range(x.shape[0]-1):
    y[i+1,:] = y[i,:] + stats.norm.rvs(mu,sigma,size=5)

for i in range(y.shape[1]):
    plt.plot(x,y[:,i])

plt.show()

・実行結果

上記の生成にあたっては「標準ブラウン運動」と同じ乱数を用いたので、$\mu=0.2$に基づく増加のトレンド以外は類似の動きであることも確認できる。

確率過程のパラメータ推定

ブラウン運動$B_t \sim N(\mu t, \sigma^2 t)$に対して、$Z_k = B_{k \Delta}-B_{(k-1) \Delta} \sim N(\mu \Delta, \sigma^2 \Delta), k=1,2,…,n$とおくとき、$\mu, \sigma^2$の推定値$\hat{\mu}, \hat{\sigma^2}$は下記のように表せる。
$$
\large
\begin{align}
\hat{\mu} \Delta &= \frac{1}{n} \sum_{k=1}^{n} Z_k \\
\hat{\sigma}^2 \Delta + (\hat{\mu} \Delta)^2 &= \frac{1}{n} \sum_{k=1}^{n} Z_k^2
\end{align}
$$

パラメータ推定にあたって、最尤法・モーメント法のどちらを用いても上記の式が導出できる。上記を$\hat{\mu}, \hat{\sigma}^2$に関して解くと下記が得られる。
$$
\large
\begin{align}
\hat{\mu} &= \frac{1}{n \Delta} \sum_{k=1}^{n} Z_k \\
\hat{\sigma}^2 &= \frac{1}{n \Delta} \sum_{k=1}^{n} Z_k^2 – \hat{\mu}^2 \Delta
\end{align}
$$

ここまでの内容を元に前項で生成を行った確率過程のパラメータ推定を下記を実行することで行うことができる。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(0)

mu, sigma = 0.2, 1
x = np.arange(0,100,1)
y = np.zeros([x.shape[0],5])

y[0,:] = stats.norm.rvs(size=5)
for i in range(x.shape[0]-1):
    y[i+1,:] = y[i,:] + stats.norm.rvs(mu,sigma,size=5)

diff = y
diff[1:diff.shape[0],:] = diff[1:diff.shape[0],:] - diff[0:(diff.shape[0]-1),:]

mu_estimate = np.mean(diff, axis=0)
sigma2_estimate = np.mean(diff**2, axis=0) - mu_estimate**2

print("Estimated mu: {}".format(mu_estimate))
print("Estimated sigma^2: {}".format(sigma2_estimate))

・実行結果

> print("Estimated mu: {}".format(mu_estimate))
Estimated mu: [ 0.09504339  0.1811264   0.13958733  0.10372601  0.34374467]
> print("Estimated sigma^2: {}".format(sigma2_estimate))
Estimated sigma^2: [ 1.01006428  0.91386789  0.99829055  1.06096954  0.92964815]

サンプル数を$n=100$から$n=10,000$に増やし、推定値の変化を確認する。

np.random.seed(0)

mu, sigma = 0.2, 1
x = np.arange(0,10000,1)
y = np.zeros([x.shape[0],5])

y[0,:] = stats.norm.rvs(size=5)
for i in range(x.shape[0]-1):
    y[i+1,:] = y[i,:] + stats.norm.rvs(mu,sigma,size=5)

diff = y
diff[1:diff.shape[0],:] = diff[1:diff.shape[0],:] - diff[0:(diff.shape[0]-1),:]

mu_estimate = np.mean(diff, axis=0)
sigma2_estimate = np.mean(diff**2, axis=0) - mu_estimate**2

print("Estimated mu: {}".format(mu_estimate))
print("Estimated sigma^2: {}".format(sigma2_estimate))

・実行結果

> print("Estimated mu: {}".format(mu_estimate))
Estimated mu: [ 0.1966156   0.2103942   0.18599239  0.1893911   0.19853666]
> print("Estimated sigma^2: {}".format(sigma2_estimate))
Estimated sigma^2: [ 1.00373888  1.00122736  0.99122914  0.99615878  0.97360955]

サンプル数を$n=10,000$に設定すると、$\mu=0.2$周辺の推定値が確認でき、この推定値は概ね妥当であると考えることができる。

また、ここでは$\Delta=1$を元に考えたが、高頻度観測(high frequent observations)を行う際も$\Delta$の値を変えるだけで、同様に下記の式を適用し、パラメータ推定を行うことができることも抑えておくと良い。
$$
\large
\begin{align}
\hat{\mu} &= \frac{1}{n \Delta} \sum_{k=1}^{n} Z_k \\
\hat{\sigma}^2 &= \frac{1}{n \Delta} \sum_{k=1}^{n} Z_k^2 – \hat{\mu}^2 \Delta
\end{align}
$$

ポアソン過程

ポアソン過程の生成

下記を実行することでポアソン過程$Po(\lambda t)$を生成することができる。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(0)

lamb = 0.2
x = np.arange(0,100,1)
y = np.zeros([x.shape[0],3])

y[0,:] = stats.norm.rvs(size=3)
for i in range(x.shape[0]-1):
    y[i+1,:] = y[i,:] + stats.poisson.rvs(lamb,size=3)

x_ = np.repeat(x,2)
y_ = np.zeros([200,3])
for i in range(1,200):
    y_[i,:] = y[(i-1)//2,:]

for i in range(y.shape[1]):
    plt.plot(x_,y_[:,i])

plt.show()

・実行結果

plot(x,y[:,i])をそのまま実行すると、階段状のグラフが表現できないので、x_y_を作成し、描画を行った。

ポアソン過程のパラメータ推定

ポアソン過程$N_t \sim Po(\lambda t)$に対して、$M_k = N_{k \Delta}-N_{(k-1) \Delta} \sim Po(\lambda \Delta), k=1,2,…,n$とおくとき、$\lambda$の推定値$\hat{\lambda}$は下記のように表せる。
$$
\large
\begin{align}
\hat{\lambda} = \frac{1}{n \Delta} \sum_{k=1}^{n} M_k = \frac{N_{n \Delta}}{n \Delta}
\end{align}
$$

ここまでの内容を元に前項で生成を行った確率過程のパラメータ推定を下記を実行することで行うことができる。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

np.random.seed(0)

lamb = 0.2
x = np.arange(0,100,1)
y = np.zeros([x.shape[0],3])

y[0,:] = stats.norm.rvs(size=3)
for i in range(x.shape[0]-1):
    y[i+1,:] = y[i,:] + stats.poisson.rvs(lamb,size=3)

diff = y
diff[1:diff.shape[0],:] = diff[1:diff.shape[0],:] - diff[0:(diff.shape[0]-1),:]

lambda_estimate = np.mean(diff, axis=0)

print("Estimated lambda: {}".format(lambda_estimate))

・実行結果

> print("Estimated lambda: {}".format(lambda_estimate))
Estimated lambda: [ 0.21764052  0.18400157  0.19978738]

サンプル数を$n=100$から$n=10,000$に増やし、推定値の変化を確認する。

np.random.seed(0)

lamb = 0.2
x = np.arange(0,10000,1)
y = np.zeros([x.shape[0],3])

y[0,:] = stats.norm.rvs(size=3)
for i in range(x.shape[0]-1):
    y[i+1,:] = y[i,:] + stats.poisson.rvs(lamb,size=3)

diff = y
diff[1:diff.shape[0],:] = diff[1:diff.shape[0],:] - diff[0:(diff.shape[0]-1),:]

lambda_estimate = np.mean(diff, axis=0)

print("Estimated lambda: {}".format(lambda_estimate))

・実行結果

> print("Estimated lambda: {}".format(lambda_estimate))
Estimated lambda: [ 0.20237641  0.19824002  0.20289787]

サンプル数を$n=10,000$に設定すると、$\lambda=0.2$周辺の推定値が確認でき、この推定値は概ね妥当であると考えることができる。

統計検定準1級問題解説 ~2019年6月実施 問9 標準ブラウン運動~

過去問題

過去問題は統計検定公式問題集が問題と解答例を公開しています。こちらを参照してください。


解答

[1] 解答

$\boxed{ \ \mathsf{18}\ }$ : ②

$B_t$は標準ブラウン運動であるから、
・$B_0=0$
・各$t\ge 0$に対して、$B_t\sim N(0,\ t)$
・任意の$0=t_0<t_1<\cdots<t_{n-1}<t_n$に対して、
  $B_{t_0},\ B_{t_1}-B_{t_0},\ \dots\ ,\ B_{t_n}-B_{t_{n-1}}$は互いに独立である。
・任意の$0 \ge t \gt t+h$に対して、
  $B_{t+h}-B_t$の分布は$B_h-B_0$の分布と同じである。($B_{t+h}-B_t \sim N(0,\ h)$)
・任意の与えられた実現値において、$B_t$は確率$1$で$t$に関して連続

以上の性質から、
 $B_k-B_{k-1}\sim N(0,\ 1)$
 $\therefore\ x_k-x_{k-1}=\sigma(B_k-B_{k-1})\sim N(0,\ \sigma^2)$
となり、互いに独立である。このとき、モーメント法による$\sigma$の推計値は、
$$\begin{eqnarray}\hat\sigma^2&=&\frac1n\sum_{k=1}^n(x_k-x_{k-1})^2-\left\{\frac1n\sum_{k=1}^n(x_k-x_{k-1})\right\}^2\\&=&\frac1n\sum_{k=1}^n(x_k-x_{k-1})^2\quad\left(\because\ \frac1n\sum_{k=1}^n(x_k-x_{k-1})\hat\mu=0\right)\end{eqnarray}$$観測データの$x_i$の増分の二乗の平均の計算結果から、$$\begin{eqnarray}\hat\sigma^2&=&\frac1{100}\sum_{k=1}^{100}(x_k-x_{k-1})^2=0.001224\\\therefore\quad\hat\sigma&=&\sqrt{0.001224}\fallingdotseq0.035\end{eqnarray}$$

[2] 解答

$\boxed{ \ \mathsf{19}\ }$ : ③

$[1]$と同じ理由から、
 $B_{\frac k{10}}-B_{\frac{k-1}{10}}=B_{\frac k{10}}-B_{\left(\frac k{10}-\frac1{10}\right)}\sim N(0,\ 1/10)$
 $\therefore\ x_{\frac k{10}}-x_{\frac{k-1}{10}}=\sigma(B_{\frac k{10}}-B_{\frac{k-1}{10}})\sim N(0,\ \sigma^2/10)$
となり、互いに独立である。このとき、モーメント法による$\sigma$の推計値は、$$\begin{eqnarray}\frac{\hat\sigma^2}{10}&=&\frac1n\sum_{k=1}^n(x_{\frac k{10}}-x_{\frac{k-1}{10}})^2\end{eqnarray}$$観測データの$x_i$の増分の二乗の平均の計算結果から、$$\begin{eqnarray}\frac{\hat\sigma^2}{10}&=&\frac1{1000}\sum_{k=1}^{1000}(x_{\frac k{10}}-x_{\frac{k-1}{10}})^2=0.000595\\\therefore\quad\hat\sigma&=&\sqrt{10\times0.000595}\fallingdotseq0.077\end{eqnarray}$$

[3] 解答

$\boxed{ \ \mathsf{20}\ }$ : ③

$x_{\frac k{10}}-x_{\frac{k-1}{10}}$と$y_{\frac k{10}}-y_{\frac{k-1}{10}}$は、$2$次元の独立同一な正規分布に従い、
$$\begin{eqnarray}E[x_{\frac k{10}}-x_{\frac{k-1}{10}}]=
E[y_{\frac k{10}}-y_{\frac{k-1}{10}}]&=&0\\V[x_{\frac k{10}}-x_{\frac{k-1}{10}}]&=&\sigma_1^2/10\\
V[y_{\frac k{10}}-y_{\frac{k-1}{10}}]&=&\sigma_2^2/10\\\mathrm{Cov}\left[x_{\frac k{10}}-x_{\frac{k-1}{10}},\ y_{\frac k{10}}-y_{\frac{k-1}{10}}\right]&=&\sigma_1\sigma_2\rho/10\end{eqnarray}$$
となる。このとき、モーメント法による$\sigma$の推計値は、
$$\begin{eqnarray}\frac{\hat\sigma_1^2}{10}&=&\frac1n\sum_{k=1}^n(x_{\frac k{10}}-x_{\frac{k-1}{10}})^2\\\frac{\hat\sigma_2^2}{10}&=&\frac1n\sum_{k=1}^n(y_{\frac k{10}}-y_{\frac{k-1}{10}})^2\\\frac{\hat\sigma_1\hat\sigma_2\hat\rho}{10}&=&\frac1n\sum_{k=1}^n(x_{\frac k{10}}-x_{\frac{k-1}{10}})(y_{\frac k{10}}-y_{\frac{k-1}{10}})\end{eqnarray}$$観測データの$x_i,y_i$の増分の二乗の平均と積和の平均の計算結果から、$$\begin{eqnarray}\frac{\hat\sigma_1^2}{10}&=&\frac1{1000}\sum_{k=1}^{1000}(x_{\frac k{10}}-x_{\frac{k-1}{10}})^2=0.000595\\\frac{\hat\sigma_2^2}{10}&=&\frac1{1000}\sum_{k=1}^{1000}(y_{\frac k{10}}-y_{\frac{k-1}{10}})^2=0.001008\\\frac{\hat\sigma_1\hat\sigma_2\hat\rho}{10}&=&\frac1{1000}\sum_{k=1}^{1000}(x_{\frac k{10}}-x_{\frac{k-1}{10}})(y_{\frac k{10}}-y_{\frac{k-1}{10}})=0.000292\\\therefore\quad\rho&=&\frac{0.000292}{\sqrt{0.000595\times0.001008}}\fallingdotseq0.377\end{eqnarray}$$


解説

ブラウン運動

各$t\in[0,\infty)$に対して、確率変数$B_t$が与えられたとき、その族$B=(B_t)_{t\ge0}$を確率過程という。$t$は時間を表すことが多く、確率過程は、時間とともに変動する偶発現象の数学モデルとして扱われる。
確率過程が以下の性質を持つとき、ブラウン運動という。
・$B_0=0$
・各$t\ge 0$に対して、$B_t\sim N(\mut,\ \sigma t)$(周辺分布)
・任意の$0=t_0<t_1<\cdots<t_{n-1}<t_n$に対して、確率変数$B_{t_0},\ B_{t_1}-B_{t_0},\ \dots\ ,\ B_{t_n}-B_{t_{n-1}}$は互いに独立である。(独立増分性)
・任意の$0 \ge t \gt t+h$に対して、
  $B_{t+h}-B_t$の分布は$B_h-B_0$の分布と同じである。($B_{t+h}-B_t \sim N(0,\ h)$)(定常増分性)
・任意の与えられた実現値において、$B_t$は確率$1$で$t$に関して連続

「統計学実践ワークブック」 演習問題etc Ch.12 「一般の分布に関する検定法」

当記事は「統計学実践ワークブック(学術図書出版社)」の読解サポートにあたってChapter.12の「一般の分布に関する検定法」に関して演習問題を中心に解説を行います。尤度比検定に加え、二項分布やポアソン分布に関する検定や適合度検定などは抑えておくと良いように思われました。

本章のまとめ

演習問題解説

例12.1

・標準正規分布を用いた検定
最尤推定値$\displaystyle \hat{\theta} = \frac{12}{30} = 0.4$を用いて、検定統計量の値は下記のように計算できる。
$$
\large
\begin{align}
\frac{\sqrt{n}(\hat{\theta}-\theta_0)}{\sqrt{\theta_0(1-\theta_0)}} &= \frac{\sqrt{30}(0.4-0.5)}{\sqrt{0.5 \times (1-0.5)}} \\
&= -1.0954…
\end{align}
$$

上記の絶対値が標準正規分布の上側$2.5$%点の$1.96$を下回ることから、帰無仮説は棄却できない。

・尤度比検定
尤度比は下記のように計算できる。
$$
\large
\begin{align}
2n \left( \hat{\theta} \log{\frac{\hat{\theta}}{\theta_0}} + (1-\hat{\theta}) \log{\frac{1-\hat{\theta}}{1-\theta_0}} \right) &= 60 \times \left( 0.4 \log{\frac{0.4}{0.5}} + 0.6 \log{\frac{0.6}{0.5}} \right) \\
&= 1.2081…
\end{align}
$$

ここで自由度$1$の$\chi^2$分布の上側$5$%点がおよそ$3.84$であることから、帰無仮説が棄却できないことがわかる。

例12.2

母比率の検定は、適合度検定においてカテゴリ数を$2$と考える場合に対応する。このとき、帰無仮説の$H_0: p=\theta_0$に基づく期待度数はそれぞれ$n \theta_0$と$n (1-\theta_0)$のように表せる。このとき、$\chi^2$統計量は下記のように表せる。
$$
\large
\begin{align}
\chi^2 &= \frac{(x_1 – n \theta_0)^2}{n \theta_0} + \frac{(n-x_1 – n(1-\theta_0))^2}{n (1-\theta_0)} \\
&= \frac{(1-\theta_0)(x_1 – n \theta_0)^2 + \theta_0(n \theta_0 – x_1)}{n \theta_0 (1-\theta_0)} \\
&= \frac{(x_1 – n \theta_0)^2}{n \theta_0 (1-\theta_0)} \\
&= \frac{n^2(x_1/n – \theta_0)^2}{n \theta_0 (1-\theta_0)} \\
&= \frac{n(\hat{\theta} – \theta_0)^2}{\theta_0 (1-\theta_0)}
\end{align}
$$

上記が$(12.1)$式で用いた統計量の$2$乗に対応するので、母比率に関する両側検定は適合度検定の特殊ケースであることが確認できる。

例12.3

独立性の仮説は$p_{ij} = \alpha_{i}\beta_{j}$で表せる。また、自由度は$d=(I-1)(J-1)$である。

詳細の考え方はワークブックが詳しいのでここでは省略。

問12.1

$[1]$
下記を実行することで$\chi^2$統計量を計算することができる。

import numpy as np

observed = np.array([7., 3., 5., 2., 7., 12., 13.])
estimate = np.repeat(np.mean(observed), observed.shape[0])

np.sum((observed-estimate)**2)/np.mean(observed)

・実行結果

> np.sum((observed-estimate)**2)/np.mean(observed)
15.142857142857142

$[2]$
$\chi^2_{\alpha=0.05}(7-1)=\chi^2_{\alpha=0.05}(6)=12.59$より、帰無仮説は棄却され。問い合わせの数が曜日によって異なると考えることができる。

問12.2

$[1]$
ポアソン分布は平均と分散が等しい分布であるが、ここで標準偏差の二乗が$1.7^2=2.89$のように平均の$2.99$に近い値を示すことが確認できる。これより、観測値がポアソン分布に従う可能性が示唆されると読み取れる。

$[2]$
ポアソン分布の確率関数を用いることで、下記のような計算を行うことができる。

import math
import numpy as np

def factorial(vec):
    res = np.zeros(vec.shape[0])
    for i in range(vec.shape[0]):
        res[i] = math.factorial(vec[i])
    return res

lamb, count = 2.99, 69.
x = np.arange(0,11,1)

p_x = lamb**x * np.e**(-lamb) / factorial(x)
p_x[10] = 1-np.sum(p_x[0:10])

print("prob: {}".format(p_x))
print("estimate: {}".format(p_x*count))

・実行結果

> print("prob: {}".format(p_x))
prob: [ 0.05028744  0.15035944  0.22478736  0.22403807  0.16746845  0.10014614
  0.04990616  0.02131706  0.00796725  0.0026469   0.00079142]
> print("estimate: {}".format(p_x*count))
estimate: [  3.46983313  10.37480107  15.5103276   15.45862651  11.55532331
   6.91008334   3.44352487   1.47087705   0.5497403    0.18263594
   0.07422687]

上記の結果を元に表を作成することができる。

$[3]$
下記の計算に基づいて$\chi^2$統計量が$16.37$であることを示すことができる。

import math
import numpy as np

def factorial(vec):
    res = np.zeros(vec.shape[0])
    for i in range(vec.shape[0]):
        res[i] = math.factorial(vec[i])
    return res

lamb, count = 2.99, 69.
x = np.arange(0,11,1)
observed = np.array([4., 7., 17., 18., 12., 7., 3., 0., 0., 0., 1.])

p_x = lamb**x * np.e**(-lamb) / factorial(x)
p_x[10] = 1-np.sum(p_x[0:10])
estimate = p_x*count
chi2_1 = np.sum((observed-estimate)**2/estimate)
estimate[10] = 0.07
chi2_2 = np.sum((observed-estimate)**2/estimate)

print("chi2-statistic 1: {}".format(chi2_1))
print("chi2-statistic 2: {}".format(chi2_2))

・実行結果

> print("chi2-statistic 1: {}".format(chi2_1))
chi2-statistic 1: 15.5647586322
> print("chi2-statistic 2: {}".format(chi2_2))
chi2-statistic 2: 16.3740364649

上記のプログラムのメインの処理では$15.56$が得られたが、estimate[10]に$0.07$を代入すると$16.37$が得られることが確認できた。これは誤差の丸め方によって変わる。

また、ここで棄却域を考えるにあたっては自由度が$11-1-1=9$より、$\chi^2_{\alpha=0.05}(9)=16.92$を用いる。統計量を大きくしている理由に関しては、下記のように$10$回以上の期待度数が小さいことに起因する。

print((1.-0.07)**2/0.07)

・実行結果

> print((1.-0.07)**2/0.07)
12.3557142857

$[4]$
下記を実行することで。$\chi^2$統計量を計算できる。

import math
import numpy as np

def factorial(vec):
    res = np.zeros(vec.shape[0])
    for i in range(vec.shape[0]):
        res[i] = math.factorial(vec[i])
    return res

lamb, count = 2.99, 69.
x = np.arange(0,7,1)
observed = np.array([4., 7., 17., 18., 12., 7., 4.])

p_x = lamb**x * np.e**(-lamb) / factorial(x)
p_x[6] = 1-np.sum(p_x[0:6])
estimate = p_x*count
chi2 = np.sum((observed-estimate)**2/estimate)

print("chi2-statistic: {}".format(chi2))

・実行結果

> print("chi2-statistic: {}".format(chi2))
chi2-statistic: 2.27565946483

上記は$\chi^2_{\alpha=0.10}(7-1-1)=9.24$よりも小さい。一方で$[3]$の結果の$16.37$は$\chi^2_{\alpha=0.10}(11-1-1)=14.68$よりも大きいので、$P$値で比較すると、上陸数をまとめたときの方があてはまりが良いと考えることができる。

問12.3

$$
\large
\begin{align}
\frac{\hat{\theta}_1-\hat{\theta}_2}{\sqrt{\left(\frac{1}{n_1}+\frac{1}{n_2}\right) \hat{\theta}_{*}(1-\hat{\theta}_{*})}} \geq z_{\alpha} \quad (12.4)
\end{align}
$$

上記で表した$(12.4)$式の左辺に対し、$n_1=114, n_2=107, \hat{\theta}_1=40/114, \hat{\theta}_2=62/107, \hat{\theta}_{*}=102/221$を代入し、計算を行う。

import numpy as np

t_statistic = (40./114.-62./107.)/np.sqrt((1./114.+1./107)*(102./221.)*(1.-102./221.))

print("t_statistic: {:.3f}".format(t_statistic))

・実行結果

> print("t_statistic: {:.3f}".format(t_statistic))
t_statistic: -3.406

参考

・準1級関連まとめ
https://www.hello-statisticians.com/toukeikentei-semi1

「統計学実践ワークブック」 演習問題etc Ch.13 「ノンパラメトリック法」

当記事は「統計学実践ワークブック(学術図書出版社)」の読解サポートにあたってChapter.13の「ノンパラメトリック法」に関して演習問題を中心に解説を行います。ノンパラメトリック法はサンプルサイズが小さい際に有効とされるケースが多く、演習を通して抑えておくと良いと思われました。

本章のまとめ

ウィルコクソンの順位和検定

「ノンパラメトリック法による仮説検定」の「ウィルコクソンの順位和検定」で取り扱った。

並べ替え検定

「ノンパラメトリック法による仮説検定」の「並べ替え検定」で取り扱った。

符号付き順位検定

「ノンパラメトリック法による仮説検定」の「符号付き順位検定」で取り扱った。

演習問題解説

例13.1

A群、B群にそれぞれ小さい順から順位を割り振り、順位和を計算すると下記のようにまとめられる。
$$
\large
\begin{array}{|c|*5{c|}}\hline \mathrm{Set} & \mathrm{Rank} & … & … & … & \mathrm{sum rank} \\
\hline A & 2 & 1 & 6 & & 9 \\
\hline B & 4 & 5 & 3 & 7 & 19 \\
\hline
\end{array}
$$

ここでA群に着目すると、選び方は${}_{7} C_{3}=35$通りある。この$35$通りのうち、順位和が$9$以下になる組み合わせは下記のように表せる。
$$
\large
\begin{align}
& 1,2,3 \qquad 1,2,4 \qquad 1,2,5 \\
& 1,2,6 \qquad 1,3,4 \qquad 1,3,5 \\
& 2,3,4
\end{align}
$$

よって、ここでの片側$P$値は$7/35=1/5=0.2$のように計算できる。

例13.2

下記を実行することで並べ替え検定における片側$P$値を得ることができる。

import numpy as np

x = np.array([30., 20., 52., 40., 50., 35., 60.])
ave_x_A = np.mean(x[0:3])

count_A = 0.
count_all = 0.
for i in range(7-2):
    for j in range(i+1,7-1):
        for k in range(j+1,7):
            count_all += 1.
            if (x[i]+x[j]+x[k])/3. <= ave_x_A:
                count_A += 1.

print("P-value: {}/{}={}".format(int(count_A), int(count_all), count_A/count_all))

・実行結果

> print("P-value: {}/{}={}".format(int(count_A), int(count_all), count_A/count_all))
P-value: 5/35=0.142857142857

並べ替え検定における片側$P$値の計算は少々大変なので、ここで用いたように何らかのプログラムを元に計算を行うのが良いと思われる。

例13.3

下記を実行することで検定を行うことができる。

import numpy as np
from scipy import stats

n, t = 35., 420.
mean_t = n*(n+1)/4.
var_t = n*(n+1)*(2*n+1)/24.

z = (t-mean_t)/np.sqrt(var_t)

if z > stats.norm.ppf(1-0.05):
    print("z: {:.2f}, P-value: {:.3f}, reject H_0.".format(z,1.-stats.norm.cdf(z)))
else:
    print("z: {:.2f}, P-value: {:.3f}, accept H_0.".format(z,1.-stats.norm.cdf(z)))

・実行結果

z: 1.72, P-value: 0.043, reject H_0.

・解説
$n(n+1)/4$は$1$から$n$までの和の$n(n+1)/2$を2で割った値であることも抑えておくと良いと思います。

例13.4

下記のような計算を実行することでクラスカル・ウォリス検定の検定統計量と、対応する$P$値の計算を行うことができる。

import numpy as np
from scipy import stats

n = np.array([5., 5., 5., 5.])
N, N_med = np.sum(n), (np.sum(n)+1.)/2.
R = np.array([15., 14., 7., 6.])

H_stat = (12.*np.sum(n*(R-N_med)**2))/(N*(N+1))

print("H-statistic: {:.2f}".format(H_stat))
print("P-value: {}".format(1.-stats.chi2.cdf(H_stat,3)))

・実行結果

> print("H-statistic: {:.2f}".format(H_stat))
H-statistic: 9.29
> print("P-value: {}".format(1.-stats.chi2.cdf(H_stat,3)))
P-value: 0.0257237409592

上記の計算結果より、有意水準$5$%で帰無仮説は棄却でき、各群に差があるといえる。

例13.5

下記を実行することで、スピアマンの順位相関係数とケンドールの順位相関係数を計算することができる。

import numpy as np

x = np.array([1., 2., 3., 4., 5., 6., 7.])
y = np.array([1., 3., 2., 6., 4., 5., 7.])

P, N = 0., 0.
for i in range(7):
    for j in range(i+1,7):
        if (x[i]-x[j])*(y[i]-y[j])>0:
            P += 1.
        else:
            N += 1.

cor_spear = 1. - 6.*np.sum((x-y)**2)/(7.*(7.**2-1))
cor_kendall = 2*(P-N)/(7.*(7.-1.))

print("Spearman correlation coef: {}".format(cor_spear))
print("Kendall correlation coef: {}".format(cor_kendall))

・実行結果

> print("Spearman correlation coef: {}".format(cor_spear))
Spearman correlation coef: 0.857142857143
> print("Kendall correlation coef: {}".format(cor_kendall))
Kendall correlation coef: 0.714285714286

問13.1

$[1]$
下記で取り扱った。
https://www.hello-statisticians.com/toukei-kentei-semi1-kakomon-cat/kakomon-semi1-2018/stat_certifi_semi1_18_5.html#2

$[2]$
下記で取り扱った。
https://www.hello-statisticians.com/toukei-kentei-semi1-kakomon-cat/kakomon-semi1-2018/stat_certifi_semi1_18_5.html#3

$[3]$
小さい順に符号付き順位を割り振ると、$-1,3,2$が得られる。これは全$2^3=8$通りのうち、$1,3,2$の次に符号付き順位和が大きくなるので、対応する$P$値は$2/8=0.25$であることがわかる。

$[4]$
$\displaystyle \frac{1}{2^n} < 0.05$を$n$に関して解けばよい。不等号は$2^n>20$と同じであるが、$n$が自然数であることを考慮すると、$n \geq 5$であることがわかる。

参考

・準1級関連まとめ
https://www.hello-statisticians.com/toukeikentei-semi1

「統計学実践ワークブック」 演習問題etc Ch.25 「因子分析・グラフィカルモデル」

当記事は「統計学実践ワークブック(学術図書出版社)」の読解サポートにあたってChapter.25の「主成分分析」に関して演習問題を中心に解説を行います。実際の分析にあたって用いられる場合がある手法なので、演習を通して抑えておくと良いと思われました。

本章のまとめ

$1$因子の因子分析(factor analysis)

因子分析(factor analysis)は$p$個の変数によって表される現象を少数の共通要因によって単純化を行う手法である。このとき、共通の要因を共通因子(common factor)と呼ぶ。

数学の試験を例に以下具体的に考える。生徒$i$の問$j$における得点を$x_{ij}$とおく。このとき$x_{ij}$を生徒$i$の因子スコア$f_i \sim N(0,1)$などを元に下記のように考えられると仮定する。
$$
\large
\begin{align}
x_{ij} = a_j f_{i} + d_j u_{ij}
\end{align}
$$
上記は$x_{ij}$が$1$つの共通因子$f_i$の関数で表されることから、$1$因子モデル(single factor model)と言われる。

ここで「共通因子$f_i$と各問$j$の関係を表す係数$a_j$」は「因子負荷量(factor loading)」と呼ばれ、共通因子では説明できない部分の$d_j u_{ij}$に関しては$d_j$を独自係数(coefficient unique factor)、$u_{ij}$を独自性因子(unique factor)と呼ぶ。

また、$x_{ij}$に関しては$V[x_{ij}]=a_j^2+d_j^2$が成立し、$a_j^2$を問$j$の共通性(communality)、$d_j^2$を問$j$の独自性(uniqueness)ということも抑えておくと良い。

実際に用いるにあたっては観測可能なのが$x_{ij}$だけであるので、$a_j, f_{i}, d_j, u_{ij}$はそれぞれ最尤推定などを用いて推定を行う。

多因子の因子分析

一般に共通因子が$K$個の場合は$K$因子モデルと言われ、下記のように表される。
$$
\large
\begin{align}
x_{ij} = a_{j1} f_{i1} + a_{j2} f_{i2} + … + a_{jK} f_{iK} + d_j u_{ij}
\end{align}
$$

多因子の因子分析を行うにあたっては、因子負荷量$a_{jk}$に関する行列の因子負荷行列の見方を抑えておくと良い。

統計学実践ワークブック表25.1より作成

上記のような推定値が得られたとき、「方程式」、「微積分」、「複素数」は問1と問2のどちらの因子負荷量も大きいが、「三角関数」と「平面図形」は問1に対応する因子負荷量がマイナスである一方で問2の因子負荷量がプラスであることが確認できる。

これにより、因子の効果を考察することができるが、このままだとわかりにくいので回転を考えることで解釈を行いやすくする手法がある。

統計学実践ワークブック表25.1より作成

たとえば上記のように変換を行えば第$1$因子で「方程式」、「微積分」、「複素数」、第$2$因子で「三角関数」と「平面図形」をそれぞれ考慮すれば良いことから解釈が行いやすい。このときに行う手続きを回転(rotation)という。

回転に関しては、主に「バリマックス法(varimax rotation)」などの「直交回転(orghogonal rotation)」や、直交しない回転である「斜交回転(oblique rotation)」などの考え方があることを抑えておくと良い。

構造方程式(Stractural equation)

統計学実践ワークブック図25.1より作成

$1$因子モデルは上記のように表すこともできるが、上記のような有向グラフに基づいて考えるモデルを構造方程式モデル(Stractural equation model)という。よって、「因子分析モデル」は「構造方程式モデル」の一種であると考えることができる。

統計学実践ワークブック図25.3より作成

より具体的に考えるにあたって、上記を元に下記のような数式を考える。
$$
\large
\begin{align}
X_{2} &= a X_{1} + u \\
X_{3} &= b X_{2} + v \\
X_{1}, &X_{2}, X_{3} \sim N(0,1)
\end{align}
$$
上記の$u, v$は誤差項であり、$u$は$X_1$と$v$は$X_1, X_2$とそれぞれ無相関であると仮定する。また、回帰係数$a, b$はそれぞれパス係数(path coefficient)と呼ばれる。

ここで$1$式目の両辺に$X_1$を掛け、期待値を取ることで下記のような導出を行うことができる。
$$
\large
\begin{align}
\rho_{12} &= Cov(X_1,X_2) \\
X_{2} &= a X_{1} + u \\
X_{1} X_{2} &= a X_{1}^2 + u X_1 \\
E[X_{1} X_{2}] &= E[a X_{1}^2 + u X_1] \\
\rho_{12} &= aE[X_{1}^2] + u E[X_1] \\
&= a \\
a &= \rho_{12}
\end{align}
$$

上記では$X_1$と$X_2$の共分散を$\rho_{12}$とおいた。また、同様に$2$式目の両辺に$X_2$を掛け、期待値を取ることで下記のような導出を行うことができる。
$$
\large
\begin{align}
\rho_{23} &= Cov(X_2,X_3) \\
X_{3} &= b X_{2} + v \\
X_{2} X_{3} &= b X_{2}^2 + b X_2 \\
E[X_{2} X_{3}] &= E[b X_{2}^2 + b X_2] \\
\rho_{23} &= bE[X_{2}^2] + u E[X_2] \\
&= b \\
b &= \rho_{23}
\end{align}
$$

このようなパス係数$a, b$の導出の方法は操作変数法(instrumental variable method)と呼ばれ、様々な構造方程式に適用することができる。また、ここで確認した操作変数法の手続きは最小二乗法と等価であることも抑えておくと良い。

演習問題解説

例25.1

$[1]$
全ての因子負荷量が正であることから、共通因子は「数学の能力の高さ」と考えることができる。

$[2]$
方程式の共通性を$a_1^2$、独自性を$d_1^2$とおくと、$a_1^2+d_1^2=1$より、それぞれ下記のように計算できる。
$$
\large
\begin{align}
a_1^2 &= 0.92^2 = 0.8464 \\
d_1^2 &= 1-a_1^2 = 1-0.8464 \\
&= 0.1536
\end{align}
$$

$[3]$
因子スコアは成績の良い順と考えられるので、下記のような対応であると考えられる。

A: -1.613
B: 1.462
C: 0.193

例25.2

$$
\large
\begin{align}
X_3 = b X_2 + v
\end{align}
$$

上記で表した$(25.2)$の第$2$式に対して両辺に$X_1$をかけて、両辺の期待値を取ることで下記のように計算できる。
$$
\large
\begin{align}
E[X_3X_1] &= E[(b X_2 + v)X_1] \\
\rho_{31} &= b E[X_1X_2] + E[v]E[X_1] = b \rho_{12} \\
b &= \frac{\rho_{31}}{\rho_{12}}
\end{align}
$$

問25.1

$[1]$
因子負荷量の絶対値$|a_{12}|$は下記のように計算できる。
$$
\large
\begin{align}
|a_{12}| &= \sqrt{0.8545-(-0.92)^2} \\
&= 0.09
\end{align}
$$

$[2]$
因子$1$のスコアが小さく、因子$2$のスコアが大きいものを選べば良いのでJが適切であることがわかる。

問25.2

$[1]$
$$
\large
\begin{align}
Y &= aX + bW + u \quad (1) \\
Z &= cX + dY + v \quad (2)
\end{align}
$$

上記の$(1)$式の両辺に$X$をかけて期待値を取ると下記のように変形を行うことができる。
$$
\large
\begin{align}
Y &= aX + bW + u \\
XY &= aX^2 + bXW + uX \\
E[XY] &= aE[X^2] + bE[XW] + E[uX] \\
\rho_{xy} &= a \cdot 1 \\
a &= \rho_{xy}
\end{align}
$$

同様に$(1)$式の両辺に$W$をかけて期待値を取ると下記のように変形を行うことができる。
$$
\large
\begin{align}
Y &= aX + bW + u \\
YW &= aXW + bW^2 + uW \\
E[YW] &= aE[XW] + bE[W^2] + E[uW] \\
\rho_{yw} &= b \cdot 1 \\
b &= \rho_{yw}
\end{align}
$$

次に$(2)$式の両辺に$X, Y$をそれぞれかけて期待値を取るとそれぞれ下記のように変形できる。
$$
\large
\begin{align}
Z &= cX + dY + v \\
XZ &= cX^2 + dXY + vX \\
E[XZ] &= cE[X^2] + dE[XY] + E[vX] \\
\rho_{xz} &= c \cdot 1 + d \rho_{xy} + 0 \\
\rho_{xz} &= c + d \rho_{xy} \quad (2.1)
\end{align}
$$
$$
\large
\begin{align}
Z &= cX + dY + v \\
YZ &= cXY + dY^2 + vY \\
E[YZ] &= cE[XY] + dE[Y^2] + vY \\
\rho_{yz} &= c \rho_{xy} + d \cdot 1 + 0 \\
\rho_{yz} &= c \rho_{xy} + d \quad (2.2)
\end{align}
$$

$(2.1)$式と$(2.2)$式に関する連立方程式を解くと、$c, d$は下記のように計算することができる。
$$
\large
\begin{align}
c &= \frac{\rho_{xz}-\rho_{xy}\rho_{yz}}{1-\rho_{xy}^2} \\
d &= \frac{\rho_{yz}-\rho_{xy}\rho_{xz}}{1-\rho_{xy}^2}
\end{align}
$$

$[2]$
$[1]$の結果より$c+ad$は下記のように計算できる。
$$
\large
\begin{align}
c + ad &= \frac{\rho_{xz}-\rho_{xy}\rho_{yz}}{1-\rho_{xy}^2} + \rho_{xy} \frac{\rho_{yz}-\rho_{xy}\rho_{xz}}{1-\rho_{xy}^2} \\
&= \frac{\rho_{xz}(1-\rho_{xy}^2)}{1-\rho_{xy}^2} \\
&= \rho_{xz}
\end{align}
$$

上記より$\rho_{xz} = c + ad$を示すことができる。

$[3]$
モラルグラフは定義より下記のように描くことができる。

参考

・準1級関連まとめ
https://www.hello-statisticians.com/toukeikentei-semi1

統計検定準1級 問題解説 ~2018年6月実施 問13 メトロポリス・ヘイスティングス法~

問題

過去問題は統計検定公式が問題と解答例を公開しています。こちらを参照してください。

解答

[1] 解答

$\boxed{ \ \mathsf{記述8}\ }$ :

目標分布の確率密度関数$\pi(y), \pi(x^{(t)})$を標準密度関数$\phi(x)$を用いて表せばよい。正規分布$N(0,1)$と$N(6,1)$の確率密度関数はそれぞれ$\phi(x), \phi(x-6)$のように表されるので、採択確率$\alpha(x^{(t)},y)$は下記のように表すことができる。
$$
\large
\begin{align}
\alpha(x^{(t)},y) = \min \left( 1 , \frac{\frac{1}{4}\phi(y)+\frac{3}{4}\phi(y-6)}{\frac{1}{4}\phi(x^{(t)})+\frac{3}{4}\phi(x^{(t)}-6)} \right)
\end{align}
$$

[2] 解答

$\boxed{ \ \mathsf{記述9}\ }$ :

ステップ幅が小さいほど山の移行が起こりにくく、どちらかの山に偏った結果になりやすい。このことからア)が$a=1$、イ)が$a=6$、ウ)が$a=0.1$であることを読み取ることができる。

[3] 解答

$\boxed{ \ \mathsf{記述10}\ }$ :

繰り返し回数の少ない段階では結果が初期値に依存することから一定程度計算が終了した後からサンプリングを行う。

解説

全体的にMCMCの基本がわかっていれば解けますが、$[1]$は問題文をよく確認しないと解答できないので$[2]$や$[3]$と比べると難しいと思われました。試験中だとこのように$[1]$が難しい場合も$[2]$や$[3]$の方が簡単な場合もあるので、全体を見て解くようにすると良いと思います。

参考

・準1級関連まとめ
https://www.hello-statisticians.com/toukeikentei-semi1

統計検定準1級 問題解説 ~2018年6月実施 問12 時系列解析~

問題

過去問題は統計検定公式が問題と解答例を公開しています。こちらを参照してください。

解答

[1] 解答

$\boxed{ \ \mathsf{23}\ }$ : ⑤

図より、$6$の倍数で自己相関係数が正かつ、それ以外では負の値を取るコレログラムを選べば良い。これにより⑤が正しいことがわかる。

[2] 解答

$\boxed{ \ \mathsf{24}\ }$ : ⑤

不規則部分の分散に季節要因が残っていることから、⑤が適切でないことがわかる。

解説

$[1]$のコレログラムの読み取りや結果の解釈に関しては出題されやすいのでワークブックなどを元に抑えておくと良いように思われました。

参考

・準1級関連まとめ
https://www.hello-statisticians.com/toukeikentei-semi1