あつまれ統計の森

Hello Statisticians!

Ch.11 「乱数の性質」の章末問題の解答例　〜基礎統計学Ⅲ 自然科学の統計学(東京大学出版会)〜

投稿日: 投稿者: lib-arts

当記事は「基礎統計学Ⅲ 自然科学の統計学(東京大学出版会)」の読解サポートにあたってChapter.11の「乱数の性質」の章末問題の解説について行います。
基本的には書籍の購入者向けの解答例・解説なので、まだ入手されていない方は下記より入手をご検討ください。また、解説はあくまでサイト運営者が独自に作成したものであり、書籍の公式ページではないことにご注意ください。

自然科学の統計学 (基礎統計学)

自然科学の統計学 (基礎統計学)

3,190円(04/26 14:55時点)
Amazon

章末の演習問題について

問題11.1の解答例

Pythonの乱数生成について確認する。
https://docs.python.org/ja/3/library/random.html
上記によると周期$2^{19937}-1$のメルセンヌツイスタ(Mersenne Twister)を用いるとされている。メルセンヌツイスタは$11.2.2$節で取り扱われているM系列を元に発展させた手法で、下記が製作者の講演資料のため参考になる。

（http://www.math.sci.hiroshima-u.ac.jp/m-mat/TEACH/ichimura-sho-koen.pdf より）

実装は下記より確認できる。
http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/MT2002/CODES/mt19937ar.c

メルセンヌツイスタ法に関しては下記で詳しく取り扱った。

メルセンヌ・ツイスタ(Mersenne Twister)法まとめ

問題11.2の解答例

下記を実行することで、乗算合同法に基づき表$11.2$と同様の結果を得ることができる。

a = 16807
M = 2**31-1
X_i = 1

for i in range(15):
    X_i = X_i*a%M
    print("{}: {}".format(i+1,X_i))

・実行結果

1: 16807
2: 282475249
3: 1622650073
4: 984943658
5: 1144108930
6: 470211272
7: 101027544
8: 1457850878
9: 1458777923
10: 2007237709
11: 823564440
12: 1115438165
13: 1784484492
14: 74243042
15: 114807987

問題11.3の解答例

例$11.3$と同様な$(p,q)=(5,3)$、$2$ビットの場合の確認を行う。

p,q = 5,3
a = [1, 1, 0, 0, 1]
X, X_10 = [], []

for i in range(2*p-len(a)):
    new_a = (a[len(a)-5]+a[len(a)-3])%2
    a.append(new_a)

for i in range(p):
    X.append(a[2*i:2*(i+1)])
    X_10.append(2*a[2*i]+a[2*i+1])

for i in range(95):
    X_new_0, X_new_1 = (X[len(X)-5][0]-X[len(X)-3][0])%2, (X[len(X)-5][1]-X[len(X)-3][1])%2
    X.append([X_new_0, X_new_1])
    X_10.append(2*X_new_0+X_new_1)

print(X_10[:15])

・実行結果

> print(X_10[:15])
[3, 0, 3, 3, 2, 0, 3, 1, 3, 1, 1, 0, 0, 2, 1]

上記は本文の表$11.4$の結果に一致することが確認できる。

次に$10,000$サンプルを発生させ、適合度検定を行う。

import numpy as np

p,q = 5,3
a = [1, 1, 0, 0, 1, 1, 1, 1, 1, 0]
X, X_10 = [], []
count_X = [0, 0, 0, 0]

for i in range(p):
    X.append(a[2*i:2*(i+1)])
    X_10.append(2*a[2*i]+a[2*i+1])

for i in range(10000):
    X_new_0, X_new_1 = (X[len(X)-5][0]-X[len(X)-3][0])%2, (X[len(X)-5][1]-X[len(X)-3][1])%2
    X.append([X_new_0, X_new_1])
    X_10.append(2*X_new_0+X_new_1)
    count_X[2*X_new_0+X_new_1] += 1

print(X_10[:15])
print(count_X)
print(np.sum((np.array(count_X)-2500.)**2/2500.))

・実行結果

> print(X_10[:15])
[3, 0, 3, 3, 2, 0, 3, 1, 3, 1, 1, 0, 0, 2, 1]
> print(count_X)
[2257, 2583, 2579, 2581]
> print(np.sum((np.array(count_X)-2500.)**2/2500.))
31.496

自由度$4-1=3$の$\chi^2$分布$\chi^2(3)$の上側$5$％点の$\chi^2(3)_{\alpha=0.05}=7.815$であるので、一様性が棄却される。一方で$n=1,000$で同様の操作を行なった場合は$31.496$に対応する値が$2.84$であり、一様性は棄却できない。

問題11.4の解答例

$$
\large
\begin{align}
x_1 &= \sqrt{-2 \log{u_1}} \cdot \cos{(2 \pi u_2)} \quad (11.6)’ \\
x_2 &= \sqrt{-2 \log{u_1}} \cdot \sin{(2 \pi u_2)} \quad (11.6)’
\end{align}
$$

上記は$(11.6)$式の文字を小文字に変換し、取りまとめたものであるが、上記を$u_1$と$u_2$に関して解くことを考える。ここで$x_1^2+x_2^2$と$\displaystyle \frac{x_2}{x_1}$を考えることでそれぞれ下記のように計算できることに着目する。
$$
\large
\begin{align}
x_1^2 + x_2^2 &= (\sqrt{-2 \log{u_1}} \cdot \cos{(2 \pi u_2)})^2 + (\sqrt{-2 \log{u_1}} \cdot \sin{(2 \pi u_2)}) \\
&= -2 \log{u_1} (\cos^{2}{(2 \pi u_2)} + \sin^{2}{(2 \pi u_2)}) \\
&= -2 \log{u_1} \\
\frac{x_2}{x_1} &= \frac{\cancel{\sqrt{-2 \log{u_1}}} \cdot \sin{(2 \pi u_2)}}{\cancel{\sqrt{-2 \log{u_1}}} \cdot \cos{(2 \pi u_2)}} \\
&= \frac{\sin{(2 \pi u_2)}}{\cos{(2 \pi u_2)}} \\
&= \tan{(2 \pi u_2)}
\end{align}
$$

上記を$u_1, u_2$に関して解くと下記が得られる。
$$
\large
\begin{align}
x_1^2 + x_2^2 &= -2 \log{u_1} \\
\log{u_1} &= -\frac{x_1^2+x_2^2}{2} \\
u_1 &= \exp \left( -\frac{x_1^2+x_2^2}{2} \right) \quad (1) \\
\frac{x_2}{x_1} &= \tan{(2 \pi u_2)} \\
2 \pi u_2 &= \tan^{-1}{ \left( \frac{x_2}{x_1} \right) } \\
u_2 &= \frac{1}{2 \pi} \tan^{-1}{ \left( \frac{x_2}{x_1} \right) } \quad (2)
\end{align}
$$

上記の$(1)$式を元に$\displaystyle \frac{\partial u_1}{\partial x_1}, \frac{\partial u_1}{\partial x_2}$は下記のように計算できる。
$$
\large
\begin{align}
\frac{\partial u_1}{\partial x_1} &= \frac{\partial}{\partial x_1} \left( -\frac{x_1^2+x_2^2}{2} \right) \\
&= \exp \left( -\frac{x_1^2+x_2^2}{2} \right) \times \frac{\partial}{\partial x_1} \left( -\frac{x_1^2+x_2^2}{2} \right) \\
&= -x_1 \exp \left( -\frac{x_1^2+x_2^2}{2} \right) = -x_1u_1 \\
\frac{\partial u_1}{\partial x_2} &= \frac{\partial}{\partial x_2} \left( -\frac{x_1^2+x_2^2}{2} \right) \\
&= -x_2 \exp \left( -\frac{x_1^2+x_2^2}{2} \right) = -x_2u_1
\end{align}
$$

また、$(2)$式に対して$\displaystyle (\tan^{-1}{y})’ = \frac{1}{1+y^2}$が成立することを元に、$\displaystyle \frac{\partial u_2}{\partial x_1}, \frac{\partial u_2}{\partial x_2}$は下記のように計算できる。
$$
\large
\begin{align}
\frac{\partial u_2}{\partial x_1} &= \frac{\partial}{\partial x_1} \left( \frac{1}{2 \pi} \tan^{-1}{ \left( \frac{x_2}{x_1} \right) } \right) \\
&= \frac{1}{2 \pi} \frac{1}{\displaystyle 1+\left(\frac{x_2}{x_1}\right)^2} \frac{\partial}{\partial x_1} \left( \frac{x_2}{x_1} \right) \\
&= \frac{1}{2 \pi} \frac{1}{\displaystyle 1+\left(\frac{x_2}{x_1}\right)^2} \left( -\frac{x_2}{x_1^2} \right) \\
&= \frac{1}{2 \pi} \frac{-x_2}{x_1^2+x_2^2} \\
&= -\frac{x_2}{2 \pi(x_1^2+x_2^2)} \\
\frac{\partial u_2}{\partial x_2} &= \frac{\partial}{\partial x_2} \left( \frac{1}{2 \pi} \tan^{-1}{ \left( \frac{x_2}{x_1} \right) } \right) \\
&= \frac{1}{2 \pi} \frac{1}{\displaystyle 1+\left(\frac{x_2}{x_1}\right)^2} \left( \frac{1}{x_1} \right) \\
&= \frac{1}{2 \pi} \frac{1}{\displaystyle 1+\left(\frac{x_2}{x_1}\right)^2} \left( \frac{x_1}{x_1^2} \right) \\
&= \frac{1}{2 \pi} \frac{x_1}{x_1^2+x_2^2} \\
&= \frac{x_1}{2 \pi(x_1^2+x_2^2)}
\end{align}
$$

よって、変数変換におけるヤコビ行列式$|\det J|$は下記のように計算できる。
$$
\large
\begin{align}
|\det J| &= \left| \begin{array}{cc} \displaystyle \frac{\partial u_1}{\partial x_1} & \displaystyle \frac{\partial u_1}{\partial x_2} \\ \displaystyle \frac{\partial u_2}{\partial x_1} & \displaystyle \frac{\partial u_2}{\partial x_2} \end{array} \right| \\
&= \left| \begin{array}{cc} \displaystyle -x_1u_1 & \displaystyle -x_2u_1 \\ \displaystyle -\frac{x_2}{2 \pi(x_1^2+x_2^2)} & \displaystyle \frac{x_1}{2 \pi(x_1^2+x_2^2)} \end{array} \right| \\
&= \left| – \frac{x_1^2u_1}{2 \pi(x_1^2+x_2^2)} – \frac{x_2^2u_1}{2 \pi(x_1^2+x_2^2)} \right| \\
&= \frac{\cancel{(x_1^2+x_2^2)}u_1}{2 \pi \cancel{(x_1^2+x_2^2)}} \\
&= \frac{u_1}{2 \pi}
\end{align}
$$

したがって、$g(x_1,x_2)$に関して下記が得られる。
$$
\large
\begin{align}
g(x_1,x_2) &= f(x_1,x_2) |\det J| \\
&= \frac{u_1}{2 \pi} = \frac{1}{2 \pi} \exp \left( -\frac{x_1^2+x_2^2}{2} \right) \\
&= \frac{1}{\sqrt{2 \pi}} \exp \left( -\frac{x_1^2}{2} \right) \times \frac{1}{\sqrt{2 \pi}} \exp \left( -\frac{x_2^2}{2} \right)
\end{align}
$$

問題11.5の解答例

下記を実行することで、$X$に$1,000$個の乱数を代入することができる。

a = 11213
M = 2**31-1
X_j = 1
X_ = []

for i in range(1000):
    X_i = 0
    for j in range(12):
        X_j = X_j*a%M
        X_i += X_j/float(M)
    X_.append(X_i-6.)

X = np.array(X_)
print(X[:15])

・実行結果

> print(X[:15])
[-1.83371108 -0.72176246 -1.12172507  0.79117246 -0.89916752 -1.45735712
  1.24864932  0.86215431  0.54254827  0.93210452  0.12854572 -0.33830579
  1.3589318  -0.59776356  0.3665337 ]

ここで得たXに対して問題$5.2$と同じ要領で適合度検定を行うことを考える。

import numpy as np

cumulative = np.array([0., np.sum(X<-3.0), np.sum(X<-2.5), np.sum(X<-2.0), np.sum(X<-1.5), np.sum(X<-1.0), np.sum(X<-0.5), np.sum(X<-0.0), np.sum(X<0.5), np.sum(X<1.0), np.sum(X<1.5), np.sum(X<2.0), np.sum(X<2.5), np.sum(X<3.0), 1000.])
observed = cumulative[1:] - cumulative[:-1]
prob = np.array([0.00135, 0.00486, 0.01654, 0.04406, 0.09185, 0.14988, 0.19146, 0.19146, 0.14988, 0.09185, 0.04406, 0.01654, 0.00486, 0.00135])
estimate = prob*1000

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

・実行結果

> print(estimate)
[-1.83371108 -0.72176246 -1.12172507  0.79117246 -0.89916752 -1.45735712
  1.24864932  0.86215431  0.54254827  0.93210452  0.12854572 -0.33830579
  1.3589318  -0.59776356  0.3665337 ]
> print((observed-estimate)**2/estimate)
[ 1.35        0.26740741  2.52307134  1.13126645  0.16137725  0.68330931
  2.40536718  0.82132874  0.02358153  0.2560969   2.71637767  0.01762999
  0.26740741  1.35      ]
> print(np.sum((observed-estimate)**2/estimate))
13.9742211708

上記より$\chi^2 < 14$であり、$\chi^2 < 14 < 22.36… = \chi^2_{\alpha=0.05}(13)$より帰無仮説は棄却できない。

まとめ

Chapter.$11$の「乱数の性質」の演習について取り扱いました。

自然科学の統計学 (基礎統計学)

自然科学の統計学 (基礎統計学)

3,190円(04/26 14:55時点)
Amazon

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

下記などの内容も参考になると思います。
https://www.hello-statisticians.com/explain-terms-cat/random_sampling1.html

「Ch.11 「乱数の性質」の章末問題の解答例　〜基礎統計学Ⅲ 自然科学の統計学(東京大学出版会)〜」への2件のフィードバック

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

コメントは受け付けていません。