2020-09-03

【統計学入門（東京大学出版会）】第6章練習問題解答

統計学

東京大学出版会から出版されている統計学入門（基礎統計学Ⅰ）について第6章の練習問題の解答を書いていきます。

本章以外の解答
6.1
6.2
6.3
6.4
- i)
- ii)
6.5
6.6
- i)
- ii)
6.7
6.8
6.9
6.10

本章以外の解答

本章以外の練習問題の解答は別の記事で公開しています。
必要に応じて参照してください。

6.1

二項分布

二項分布の期待値 $E(X)$ は、

$E(X) = np$

で与えられます。

一方 $E(X^{2})$ は、

$\begin{align} E(X^{2}) &= \sum_{x} x^{2} \cdot {}_{n}C_{x} p^{x} (1-p)^{n-x} \\ &= \sum_{x} ( x ( x - 1 ) + x ) \cdot {}_{n}C_{x} p^{x} (1-p)^{n-x} \\ &= \sum_{x} x ( x - 1 ) \cdot {}_{n}C_{x} p^{x} (1-p)^{n-x} + \sum_{x} x \cdot {}_{n}C_{x} p^{x} (1-p)^{n-x} \\ &= \sum_{x} x ( x - 1 ) \cdot \frac{n(n-1)}{x(x-1)} \cdot p^{2} \cdot {}_{n-2}C_{x-2} p^{x-2} (1-p)^{n-x} + \sum_{x} x \cdot \frac{n}{x} \cdot p \cdot {}_{n-1}C_{x-1} p^{x-1} (1-p)^{n-x} \\ &= np^{2}(n-1) \cdot \sum_{x} \cdot {}_{n-2}C_{x-2} p^{x-2} (1-p)^{n-x} + np \cdot \sum_{x} {}_{n-1}C_{x-1} p^{x-1} (1-p)^{n-x} \\ &= np^{2}(n-1) \cdot 1 + np \cdot 1 \\ &= n^{2}p^{2} + np (1-p) \end{align}$

となるため、分散 $V(X)$ は、

$\begin{align} V(X) &= E(X^{2}) - E(X)^{2} \\ &= n^{2}p^{2} + np (1-p) - (np)^{2} \\ &= np(1-p) \end{align}$

となります。

ポアソン分布

ポアソン分布の期待値 $E(X)$ は、

$E(X) = \lambda$

で与えられます。

一方 $E(X^{2})$ は、

$\begin{align} E(X^{2}) &= \sum_{x} x^{2} \cdot e^{- \lambda} \cdot \frac{\lambda^{x}}{x!} \\ &= \sum_{x} ( x ( x - 1 ) + x) \cdot e^{- \lambda} \cdot\frac{\lambda^{x}}{x!} \\ &= \sum_{x} x ( x - 1 ) \cdot e^{- \lambda} \cdot \frac{\lambda^{x}}{x!} + \sum_{x} x \cdot e^{- \lambda} \cdot \frac{\lambda^{x}}{x!} \\ &= \lambda^{2} e^{- \lambda} \cdot \sum_{x} \frac{\lambda^{x 2}}{(x-2)!} + \lambda e^{- \lambda} \cdot \sum_{x} \frac{\lambda^{x-1}}{(x-1)!} \\ &= \lambda^{2} e^{- \lambda} e^{\lambda} + \lambda e^{- \lambda} e^{\lambda} = \lambda^{2} + \lambda \end{align}$

となるため、分散 $V(X)$ は、

$\begin{align} V(X) &= E(X^{2}) - E(X)^{2} \\ &= \lambda^{2} + \lambda - \lambda^{2} \\ &= \lambda \end{align}$

となります。

6.2

ポアソン分布 $P_{o}(\lambda)$ は、次の式で与えられます。

$\displaystyle P\_{o}(\lambda) = e^{- \lambda} \cdot \frac{\lambda^{x}}{x!}$

4床の空きベッドが確保されているため、ベッドが不足する確率は救急患者数が5人以上である確率を求めればよいことになります。
したがって、

$\displaystyle 1- \sum_{x=0}^{4} P\_{o}(2.5) = 1 - \sum_{x=0}^{4} e^{- 2.5} \cdot \frac{2.5^{x}}{x!}$

を求めることで答えが得られます。

上記の計算を行うPythonプログラムを次に示します。

from math import exp, pow, factorial

ans = 1.0
for x in range(5):
    ans -= exp(-2.5) * pow(2.5, x) / factorial(x)
print(ans)

上記のプログラムを実行すると、次の結果が得られます。

0.10882198108584873

6.3

負の二項分布とは、 $k$ 回目の成功を得るまでの試行回数 $N$ に関する確率分布 です。したがって最後の試行が成功となり、それ以外の $N - 1$ 回の試行では、 $k-1$ 回の成功と $x = N - k$ 回の失敗となる確率を求めればよいことになります。

成功の確率を $p$ 失敗の確率を $q$ とすると、確率分布 $f(x)$ は、

$\begin{align} f(x) &= p \cdot {}_{N-1}C_{x} p^{k-1} q^{x} \\ &= {}_{x+k-1}C_{x}p^{k} q^{x} \end{align}$

以上により、負の二項分布を導出できました。

6.4

i)

$n$ 個のコインのうち、1個のコインが表になり $n-1$ 個のコインが裏になる確率と、 $n-1$ 個のコインが表になり1個のコインが裏になる確率の和が $P$ になります。

したがって、

$\begin{align} P &= {}_{n}C_{1} p q^{n-1} + {}_{n}C_{1} p^{n-1} q \\ & = n(pq^{n-1} + p^{n-1}q) \end{align}$

ii)

繰り返し数を $x$ とすると、 $x$ 回目でi)を満たす確率 $Q$ は、

$Q = P(1-P)^{x-1}$

となるため、 $x$ の期待値 $E(X)$ は、

$\begin{align} E(X) &= \sum_{x} x \cdot Q \\ &= \sum_{x} x \cdot P(1-P)^{x-1} \end{align}$

から求めることができます。

ここで $x$ が非常に大きい（＝無限大）のときは、

$\begin{align} \sum_{x} x p q^{x-1} &= \sum_{x} ( x + 1 ) p q^{x} \\ &= \sum_{x} x p q^{x} + \sum_{x} p q^{x} \\ &= q \cdot \sum_{x} x p q^{x-1} + 1 \end{align}$

が成り立つため、

$(1-q) \cdot \sum_{x} x p q^{x-1} = 1$

の関係式が得られます。

この関係式を利用すると、

$\begin{align} E(X) &= \sum_{x} x \cdot P(1-P)^{x-1} \\ &= \frac{1}{1-(1-P)} \\ &= \frac{1}{P} \\ &= \frac{1}{n(pq^{n-1} + p^{n-1}q)} \end{align}$

が得られます。

6.5

定数 $c$

$\begin{equation} f(x) \begin{cases} \displaystyle c(1-x^{2}) & (-1 \le x \le 1) \\ 0 & (x \lt -1, 1 \lt x) \end{cases} \end{equation}$

が確率密度関数となるためには、

$\displaystyle \int_{-\infty}^{\infty} f(x) dx = 1$

を満たせばよいことになります。したがって、

$\begin{align} \int_{-\infty}^{\infty} f(x) dx &= \int_{-1}^{1} c (1 - x^{2}) dx \\ &= \frac{4}{3} c \\ &= 1 \end{align}$

より（偶関数の性質を利用）、 $c = \frac{3}{4}$ が求まります。

以降の計算では、この $c$ の値を利用して期待値などの値を求めます。すなわち、

$\begin{equation} f(x) \begin{cases} \displaystyle \frac{3}{4}(1-x^{2}) & (-1 \le x \le 1) \\ 0 & (x \lt -1, 1 \lt x) \end{cases} \end{equation}$

です。

期待値

$f(x)$ の期待値 $E(X)$ は、

$\begin{align} E(X) &= \int_{-\infty}^{\infty} x \cdot f(x) dx \\ &= \int_{-1}^{1} \frac{3}{4} x (1 - x^{2}) dx \\ & = 0 \end{align}$

となります（奇関数の性質を利用）。

分散

$\begin{align} E(X^{2}) &= \int_{-\infty}^{\infty} x^{2} \cdot f(x) dx \\ &= \int_{-1}^{1} \frac{3}{4} x^{2} (1 - x^{2}) dx \\ &= 2 \cdot \int_{0}^{1} \frac{3}{4} x^{2} (1 - x^{2}) dx \\ & = \frac{1}{5} \end{align}$

となるため、分散 $V(X)$

$\begin{align} V(X) &= E(X^{2}) - E(X)^{2} \\ &= \frac{1}{5} \end{align}$

が得られます。

歪度

$\mu = E(X) = 0$ 、 $\sigma = \sqrt{V(X)} = \frac{1}{\sqrt{5}}$ と、

$\begin{align} E(X^{3}) &= \int_{-\infty}^{\infty} x^{3} \cdot f(x) dx \\ &= \int_{-1}^{1} \frac{3}{4} x^{3} (1 - x^{2}) dx \\ & = 0 \end{align}$

より、歪度 $\alpha_{3}$ は、

$\begin{align} \alpha_{3} &= \frac{E(X^{3}) - 3 \mu E(X^{2}) + 2 \mu^{3}}{\sigma^{3}} \\ &= 0 \end{align}$

となります。

尖度

$\begin{align} E(X^{4}) &= \int_{-\infty}^{\infty} x^{4} \cdot f(x) dx \\ &= \int_{-1}^{1} \frac{3}{4} x^{4} (1 - x^{2}) dx \\ &= 2 \cdot \int_{0}^{1} \frac{3}{4} x^{4} (1 - x^{2}) dx \\ & = \frac{3}{35} \end{align}$

より、尖度 $\alpha_{4} - 3$ は、

$\begin{align} \alpha_{4} - 3 &= \frac{E(X^{4}) - 4 \mu E(X^{2}) + 6 \mu^{2} E(X^{2}) - 3 \mu^{4}}{\sigma^{4}} \\ &= \frac{\frac{3}{35}}{\frac{1}{\sqrt{5}^{4}}} -3 \\ &= \frac{6}{7} \end{align}$

となります。

6.6

i)

指数分布の確率密度関数 $f(x)$ は、次の式で与えられます（ $\lambda$ は正の値）。

$\begin{equation} f(x) \begin{cases} \lambda e^{- \lambda x} & (x \ge 0) \\ 0 & (x \lt 0) \end{cases} \end{equation}$

これを用いて、

$\begin{align} P(X \gt a + b | X \gt a) &= \frac{P(X \gt a + b \land X > a)}{P(X \gt a)} \\ &= \frac{P(X \gt a + b)}{P(X \gt a)} \\ &= \frac{\int_{a+b}^{\infty} f(x) dx}{\int_{a}^{\infty} f(x) dx} \\ &= \frac{\int_{a+b}^{\infty} \lambda e^{- \lambda x} dx}{\int_{a}^{\infty} \lambda e^{- \lambda x} dx} \\ &= \frac{e^{- \lambda (a+b)}}{e^{-\lambda a}} \\ &= e^{-\lambda b} \\ &= P(X \gt b) \end{align}$

となります。

$P(X \gt a + b | X \gt a)$ は、過去に $a$ だけの時間が過ぎた状態という前提条件をもとにして、 $b$ だけ時間を進めたときの確率を示しています。一方で $P(X \gt b)$ は、いかなる前提条件をもとにせず、 $b$ だけ時間を進めたときの確率を示しています。これらが同じ確率になっているということは、過去の時間経過がその後の確率に影響を与えていない、ということを示していると言えます。

ii)

累積分布関数 $F(x)$ は、

$\begin{align} F(x) &= \int_{-\infty}^{x} f(x) dx \\ &= \int_{0}^{x} \lambda e^{-\lambda x} dx \\ &= 1 - e^{- \lambda x} \end{align}$

となるため、

$\begin{align} \lambda (x) &= \frac{f(x)}{1-F(x)} \\ &= \frac{\lambda e^{- \lambda x}}{1 - (1 - e^{- \lambda x})} \\ &= \lambda \end{align}$

が得られます。

6.7

付表の正規分布表を利用します。

付表は上側の確率の値を示しているため、 $P(|Z| \gt c)$ の場合は、表の値の1/2となる値を見る必要があることに注意が必要です。例えば、 $P(|Z| \gt c) = 0.01$ の場合は、0.005に対応する $c$ の値を参照するといった具合です。

また本来は、内挿を考慮して値を求める必要がありますが、簡単のため2点間で近い方の値を $c$ の値として採用しています。

$P(\|Z\| \gt c)$	$c$
0.01	2.58
0.02	2.32
0.05	1.96
0.10	1.65

および

$P(Z \gt c)$	$c$
0.01	2.32
0.02	2.05
0.05	1.65
0.10	1.28

が得られます。

6.8

ベータ分布の確率密度関数 $f(x)$ は、

$\begin{equation} f(x) \begin{cases} \frac{x^{\alpha -1} (1 - x)^{\beta -1}}{B(\alpha, \beta)} & (0 \ge x \ge 1) \\ 0 & (x \lt 0, 1 \lt x) \end{cases} \end{equation}$

かつ凹関数であることから、 $f(x)$ を微分して0となる $x$ の値がモード（最頻）となります。

したがって、

$\begin{align} \frac{df(x)}{dx} &= \frac{1}{B(\alpha, \beta)} \left( (\alpha - 1) x^{\alpha - 2} (1 - x)^{\beta - 1} - (\beta - 1) x^{\alpha - 1} (1 - x)^{\beta - 2} \right) \\ &= 0 \end{align}$

を満たす $x$ を求めればよいことになります。 $B(\alpha, \beta)$ は $x$ に依存しないことに注意して計算すると、

$\displaystyle x = \frac{\alpha - 1}{\alpha + \beta - 2}$

が得られます。

なお、 $\alpha = \beta = 1$ のときはベータ分布が一様分布になることから、モードは $0 \lt x \lt 1$ の範囲で任意の値を取れる点に注意してください。

6.9

ワイブル分布の密度関数 $f(x)$ を次に示します。

$\begin{equation} f(x) \begin{cases} \displaystyle \frac{bx^{b-1}}{a^{b}} \exp \left( - \left( \frac{x}{a} \right)^{b} \right) & (x \gt 0) \\ 0 & (x \lt 0) \end{cases} \end{equation}$

累積分布関数 $P(X \le x)$ は、

$\begin{align} P(X \le x) &= \int_{-\infty}^{x} f(y) dy \\ &= \int_{0}^{x} \frac{by^{b-1}}{a^{b}} \exp \left( - \left( \frac{y}{a} \right)^{b} \right) dy \\ &= \frac{b}{a^{b}} \left[ - \frac{a^{b}}{b} \exp \left( - \frac{y^{b}}{a^{b}} \right) \right]_{y=0}^{y=x} \\ &= 1 - \exp \left( - \left( \frac{x}{a} \right)^{b} \right) \end{align}$

と求まります。ここで求めた累積分布関数は、 $x \ge 0$ を満たす場合に限定しています。 $x \lt 0$ の場合は $f(x) = 0$ となるので、累積分布関数も0になります。

6.10

標準正規分布

標準正規分布の確率密度関数 $f(x)$ は、次の式で与えられます。

$\displaystyle f(x) = \frac{1}{\sqrt{2\pi}} \exp \left( - \frac{x^{2}}{2} \right)$

したがってモーメント母関数 $M_{x}(t)$ は、変数変換 $x - t \equiv y$ とガウス積分の公式を使って求めることができます。

$\begin{align} M_{x}(t) &= \int_{-\infty}^{\infty} e^{tx} f(x) dx \\ &= \int_{-\infty}^{\infty} e^{tx} \cdot \frac{1}{\sqrt{2\pi}} \exp \left( - \frac{x^{2}}{2} \right) dx \\ &= \frac{1}{\sqrt{2\pi}} \exp \left( \frac{t^{2}}{2} \right) \int_{- \infty}^{\infty} \exp \left( - \frac{1}{2} (x - t)^{2} \right) dx \\ &= \frac{1}{\sqrt{2\pi}} \exp \left( \frac{t^{2}}{2} \right) \int_{- \infty}^{\infty} \exp \left( - \frac{y^{2}}{2} \right) dy \\ &= \frac{1}{\sqrt{2\pi}} \exp \left( \frac{t^{2}}{2} \right) \sqrt{\frac{\pi}{\frac{1}{2}}} \\ &= \exp \left( \frac{t^{2}}{2} \right) \end{align}$

ここでマクローリン展開すると、

$\begin{align} M_{x}(t) &= \exp \left( \frac{t^{2}}{2} \right) \\ &\equiv e^{y} \\ &= \sum_{n=0}^{\infty} \frac{(y - 0)^{n}}{n!} \frac{d^{n}}{dy^{n}} e^{y} |_{y=0} \\ &= \sum_{n=0}^{\infty} \frac{y^{n}}{n!} \\ &= 1 + y + \frac{y^{2}}{2!} + \frac{y^{3}}{3!} + ... \\ &= 1 + \frac{1}{2}t^{2} + \frac{1}{8}t^{4} + ... \end{align}$

となります。

一方、モーメント母関数 $M_{x}(t)$ は、

$\displaystyle M_{x}(t) = 1 + t \cdot E(X) + \frac{t^{2}}{2!} \cdot E(X^{2}) + \frac{t^{3}}{3!} \cdot E(X^{3}) + \frac{t^{4}}{4!} \cdot E(X^{4}) + ...$

という性質があるため、

$E(X) = \mu = 0$

$E(X^{2}) = 1$

$E(X^{3}) = 0$

$E(X^{4}) = 3$

$V(X) = \sigma^{2} = E(X^{2}) - E(X)^{2} = 1$

が得られます。

よって尖度 $\alpha_{4} - 3$ は、

$\begin{align} \alpha_{4} - 3 &= \frac{E(X^{4}) - 4 \mu E(X^{3}) + 6 \mu^{2} E(X^{2}) - 3 \mu^{4}}{\sigma^{4}} - 3 \\ &= \frac{3}{1} - 3 \\ &= 0 \end{align}$

となります。

指数分布

指数分布の確率密度関数 $f(x)$ は、次の式で与えられます。

$\begin{equation} f(x) \begin{cases} \displaystyle \lambda e^{- \lambda x} & (x \ge 0) \\ 0 & (x \lt 0) \end{cases} \end{equation}$

したがってモーメント母関数 $M_{x}(t)$ は、次のようになります。なお、 $\lambda > t$ とします。

$\begin{align} M_{x}(t) &= \int_{-\infty}^{\infty} e^{tx} f(x) dx \\ &= \int_{0}^{\infty} e^{tx} \cdot \lambda e^{- \lambda x} dx \\ &= \frac{\lambda}{\lambda -t} \end{align}$

ここでマクローリン展開すると、

$\begin{align} M_{x}(t) &= \frac{\lambda}{\lambda -t} \\ &= \sum_{n=0}^{\infty} \frac{(t - 0)^{n}}{n!} \frac{d^{n}}{dt^{n}} \frac{\lambda}{\lambda -t} |_{t=0} \\ &= 1 + \frac{1}{\lambda}t + \frac{1}{\lambda^{2}}t^{2} + \frac{1}{\lambda^{3}}t^{3} + \frac{1}{\lambda^{4}}t^{4} + ... \end{align}$

となります。

したがって、

$\displaystyle E(X) = \mu = \frac{1}{\lambda}$

$\displaystyle E(X^{2}) = \frac{2!}{\lambda}$

$\displaystyle E(X^{3}) = \frac{3!}{\lambda}$

$\displaystyle E(X^{4}) = \frac{4!}{\lambda}$

$\displaystyle V(X) = \sigma^{2} = E(X^{2}) - E(X)^{2} = \frac{1}{\lambda^{2}}$

が得られます。

よって尖度 $\alpha_{4} - 3$ は、

$\begin{align} \alpha_{4} - 3 &= \frac{E(X^{4}) - 4 \mu E(X^{3}) + 6 \mu^{2} E(X^{2}) - 3 \mu^{4}}{\sigma^{4}} - 3 \\ &= 9 - 3 \\ &= 6 \end{align}$

となります。

2020-09-01

【Python】機械学習ライブラリにおける乱数固定

Python 機械学習

機械学習向けに提供されるライブラリは、内部で乱数を使っていることが多く、何も考えずにプログラムを書くとプログラム実行の度に結果が変わってしまいます（再現性がありません）。再現性のあるプログラムを書くためには乱数を固定化する必要がありますが、利用するライブラリによって乱数を固定化する方法が異なります。

本記事では、自分が日頃利用しているライブラリについて、再現性を確保するために必要となる手続きをまとめたいと思います。なお、本記事は自分が頻繁に利用するライブラリに絞って説明していますので、すべての機械学習のライブラリについて網羅できているわけではありません。また、本記事は随時更新していく予定です。

Python

ハッシュ生成のランダム化を無効化参考

os.environ["PYTHONHASHSEED"] = str(seed)

random

乱数Seed固定参考

random.seed(seed)

numpy

乱数Seed固定参考

np.random.seed(seed)

PyTorch

乱数Seed固定参考

torch.manual_seed(seed)

CUDA専用の乱数Seed固定API torch.cuda.manual_seed もありますが、上記のAPIを実行するだけでCUDA側の乱数Seedも固定してくれるようです。

cudnn内の非決定的な処理の固定化参考

torch.backends.cudnn.deterministic = True
torch.backends.cudnn.benchmark = False

torch.backends.cudnn.benchmark に False にすると最適化による実行の高速化の恩恵は得られませんが、テストやデバッグ等に費やす時間を考えると結果としてトータルの時間は節約できる、と公式のドキュメントには記載されていました。

2020-08-31

【統計学入門（東京大学出版会）】第5章練習問題解答

統計学

東京大学出版会から出版されている統計学入門（基礎統計学Ⅰ）について第5章の練習問題の解答を書いていきます。

本章以外の解答
5.1
- i)
- ii)
- iii)
5.2
5.3
- i)
- ii)
5.4
5.5
5.6
5.7
5.8

本章以外の解答

本章以外の練習問題の解答は別の記事で公開しています。
必要に応じて参照してください。

5.1

i)

密度関数

一様分布の密度関数 $f(x)$ を次に示します。

$\begin{equation} f(x) \begin{cases} \displaystyle \frac{1}{6} & (0 \le x \le 6) \\ 0 & (x \lt 0, 6 \lt x) \end{cases} \end{equation}$

期待値

期待値 $E(X)$ は、

$\begin{align} \displaystyle E(X) &= \int_{-\infty}^{\infty} x \cdot f(x) dx \\ &= \int_{0}^{6} \frac{x}{6} dx \\ &= 3 \end{align}$

となります。

分散

$\begin{align} \displaystyle E(X^{2}) &= \int_{-\infty}^{\infty} x^{2} \cdot f(x) dx \\ &= \int_{0}^{6} \frac{x^{2}}{6} dx \\ &= 12 \end{align}$

と求まるので、

分散は $V(X)$ は、

$\begin{align} \displaystyle V(X) &= E(X^{2}) - E(X)^{2} \\ &= 12 - 3^{2} \\ &= 3 \end{align}$

となります。

ii)

チェビシェフの不等式は、次の式で与えられます。

$\displaystyle P(|X - \mu| \ge k \sigma) \le \frac{1}{k^{2}}$

したがって、i)の結果を利用して、

$\displaystyle P(|X - 3| \ge \sqrt{3} k) = P(3 - \sqrt{3} k \ge X \ge 3 + \sqrt{3} k) \le \frac{1}{k^{2}}$

が成り立つことを適当なk（ここでは $k = 1$ とします）について確認すればよいことになります。

$\begin{align} \displaystyle P(3 - \sqrt{3} k \ge X \ge 3 + \sqrt{3} k) &= P(3 - \sqrt{3} \ge X \ge 3 + \sqrt{3}) \\ &= \int_{3 - \sqrt{3}}^{3 + \sqrt{3}} f(x) dx \\ &= \int_{3 - \sqrt{3}}^{3 + \sqrt{3}} \frac{1}{6} dx \\ &= 0.57735 \quad (\le 1) \end{align}$

となり、不等式が成り立つことが確認できました。

iii)

密度関数 $f(x)$ を次に示します。

$\begin{equation} f(x) \begin{cases} \displaystyle 1 & (0 \le x \le 1) \\ 0 & (x \lt 0, 1 \lt x) \end{cases} \end{equation}$

これから、

$\begin{align} \displaystyle E(X) &= \int_{-\infty}^{\infty} x \cdot f(x) dx \\ &= \int_{0}^{1} x dx \\ &= \frac{1}{2} \end{align}$

$\begin{align} \displaystyle E(X^{2}) &= \int_{-\infty}^{\infty} x^{2} \cdot f(x) dx \\ &= \int_{0}^{1} x^{2} dx \\ &= \frac{1}{3} \end{align}$

$\begin{align} \displaystyle E(X^{3}) &= \int_{-\infty}^{\infty} x^{3} \cdot f(x) dx \\ &= \int_{0}^{1} x^{3} dx \\ &= \frac{1}{4} \end{align}$

$\begin{align} \displaystyle E(X^{4}) &= \int_{-\infty}^{\infty} x^{4} \cdot f(x) dx \\ &= \int_{0}^{1} x^{4} dx \\ &= \frac{1}{5} \end{align}$

と求まります。

歪度

歪度 $\alpha_{3}$ は、

$\displaystyle \alpha_{3} = \frac{E((X - \mu)^{3})}{\sigma^{3}}$

で求めますが、

$\begin{align} \displaystyle E((X - \mu)^{3}) &= E(X^{3}) - 3 \mu E(X^{2}) + 2 \mu^{2} \\ &= \frac{1}{4} - 3 \cdot \frac{1}{2} \cdot \frac{1}{3} + 2 \cdot \left( \frac{1}{2} \right)^{3} \\ &= 0 \end{align}$

であるため、 $\alpha_{3} = 0$ となります。

尖度

尖度 $\alpha_{4} -3$ は、

$\displaystyle \alpha_{4} - 3 = \frac{E((X - \mu)^{4})}{\sigma^{4}} - 3$

で求めますが、

$\begin{align} \displaystyle \sigma^{2} &= V(X) \\ &= E(X^{2}) - E(X)^{2} \\ &= \frac{1}{3} - \left( \frac{1}{2} \right)^{2} \\ &= \frac{1}{12} \end{align}$

$\begin{align} \displaystyle E((X - \mu)^{4}) &= E(X^{4}) - 4 \mu E(X^{3}) + 6 \mu^{2} E(X^{2}) - 3 \mu^{4} \\ &= \frac{1}{5} - 4 \cdot \frac{1}{2} \cdot \frac{1}{4} +6 \cdot \left( \frac{1}{2} \right)^{2} \cdot \frac{1}{3} - 3 \cdot \left( \frac{1}{2} \right)^{4} \\ &= \frac{1}{80} \end{align}$

から、

$\begin{align} \displaystyle \alpha_{4} - 3 &= \frac{E((X - \mu)^{4})}{\sigma^{4}} - 3 \\ &= \frac{\frac{1}{80}}{\left( \frac{1}{12} \right)^{2}} - 3 \\ &= - \frac{6}{5} \end{align}$

となります。

5.2

各等級の当選金を $x_{i}$ 、本数を $n_{i}$ とすると、宝くじの期待値 $E(X)$ は、次の式で表されます。
なお、全本数を $N$ としました。

$\displaystyle E(X) = \sum_{i} x_{i} \cdot \frac{n_{i}}{N}$

上記の式を用いて期待値を求めるPythonプログラムを次に示します。

import numpy as np

money = np.array([
    40000000,
    10000000,
    200000,
    10000000,
    100000,
    1000000,
    140000,
    10000,
    1000,
    200
])
num = np.array([
    7,
    14,
    903,
    5,
    645,
    130,
    130,
    1300,
    26000,
    1300000
])

total = np.sum(money * num) / 13000000
print(total)

上記のプログラムを実行した結果を次に示します。

89.4076923076923

5.3

i)

n回目に初めてコインが表になる確率 $P(X)$ は、n-1回目まで裏が出ることと同じであるため、次の式で表されます。

$\begin{align} \displaystyle P(X) &= \left( \frac{1}{2} \right)^{n-1} \cdot \frac{1}{2} \\ &= \frac{1}{2^{n}} \quad (n = 1, 2, 3, ...) \end{align}$

ii)

i)の結果を用いて期待値 $E(X)$ を求めると、

$\begin{align} \displaystyle E(X) &= \sum_{n=1}^{\infty} x_{n} P(X=x_{n}) \\ &= \sum_{n=1}^{\infty} 2^{n} \cdot \frac{1}{2^{n}} \\ &= \sum_{n=1}^{\infty} 1 \\ &= \infty \end{align}$

となり題意は示されました。

5.4

$E(X-a)^{2}$ を式変形すると、

$E(X-a)^{2} = E(X^{2}) - 2aE(X) + a^{2}$

となります。上記の式を最小にするためには、 $E(X-a)^{2}$ を $a$ について微分した値が0になればよいので、

$\begin{align} \displaystyle \frac{dE}{da} &= -2E(X) + 2a \\ &= 0 \end{align}$

したがって、 $a = E(X)$ のとき $E(X-a)^{2}$ が最小値 $E(X^{2}) - E(X)^{2}$ をとります。

5.5

正n面体で各面が出る確率は、 $\frac{1}{n}$ です。ここから、各期待値と分散が求まります。

期待値

期待値 $E(X)$ は、1次式の和の公式 $\sum_{i=1}^{n} = \frac{1}{2} n(n+1)$ を用いて、

$\begin{align} \displaystyle E(X) &= \sum_{i=1}^{n} x_{i} \cdot \frac{1}{n} \\ &= \frac{1}{2} n(n+1) \cdot \frac{1}{n} \\ &= \frac{n+1}{2} \end{align}$

となります。

分散

分散 $V(X)$ は、1次式の和の公式と2次式の和の公式 $\sum_{i=1}^{n} = \frac{1}{6} n(n+1)(2n+1)$ を用いて、

$\begin{align} \displaystyle E(X^{2}) &= \sum_{i=1}^{n} x_{i}^{2} \cdot \frac{1}{n} \\ &=\frac{1}{6} n(n+1)(2n+1) \cdot \frac{1}{n} \\ &= \frac{(n+1)(2n+1)}{6} \end{align}$

$\begin{align} \displaystyle V(X) &= E(X^{2}) - E(X)^{2} \\ &= \frac{(n+1)(2n+1)}{6} - \left( \frac{n+1}{2} \right)^{2} \\ &= \frac{n^{2} - 1}{12} \end{align}$

となります。

5.6

確率変数 $X$ の密度関数 $f(x)$ は、次のように与えられます。

$\begin{equation} f(x) \begin{cases} 1 & (0 \le x \le 1) \\ 0 & (x \lt 0, 1 \lt x) \end{cases} \end{equation}$

累積分布関数

$X^{2}$ の累積分布関数 $P(X^{2} \le y)$ は、

$\begin{align} \displaystyle P(X^{2} \le y) &= P(- \sqrt{y} \le X \le \sqrt{y}) \\ &= \int_{-\sqrt{y}}^{\sqrt{y}} f(x) dx \\ &= \int_{0}^{\sqrt{y}} 1 \cdot dx \\ &= \sqrt{y} \end{align}$

と求まります。

密度関数

密度関数 $g(y)$ は、累積分布関数を微分したものであるため、

$\begin{align} \displaystyle g(y) &= \frac{dP}{dy} \\ &= \frac{d}{dy} \sqrt{y} \\ &= \frac{1}{2\sqrt{y}} \quad (0 \le y \le 1) \end{align}$

となります。

期待値

期待値 $E(Y)$ は、

$\begin{align} \displaystyle E(Y) &= \int_{-\infty}^{\infty} y \cdot g(y) dy \\ &= \int_{0}^{1} \frac{\sqrt{y}}{2} dy \\ &= \left[ \frac{y^{3/2}}{3} \right]_{0}^{1} \\ &= \frac{1}{3} \end{align}$

と求まります。

分散

$\begin{align} \displaystyle E(Y^{2}) &= \int_{-\infty}^{\infty} y^{2} \cdot g(y) dy \\ &= \int_{0}^{1} \frac{y^{3/2}}{2} dy \\ &= \left[ \frac{y^{5/2}}{5} \right]_{0}^{1} \\ &= \frac{1}{5} \end{align}$

となるため、分散 $V(Y)$ は、

$\begin{align} \displaystyle V(Y) &= E(Y^{2}) - E(Y)^{2} \\ &= \frac{1}{5} - \left( \frac{1}{3} \right) \\ &= \frac{4}{45} \end{align}$

と求まります。

5.7

確率変数 $X$ の密度関数 $f(x)$ は、次のように与えられます。

$\displaystyle f(x) = \frac{1}{\sqrt{2\pi}} \exp \left( - \frac{x^{2}}{2} \right)$

累積分布関数

$X^{2}$ の累積分布関数 $P(X^{2} \le y)$ は、

$\begin{align} \displaystyle P(X^{2} \le y) &= P(- \sqrt{y} \le X \le \sqrt{y}) \\ &= \int_{-\sqrt{y}}^{\sqrt{y}} f(x) dx \\ &= \int_{-\sqrt{y}}^{\sqrt{y}} \frac{1}{\sqrt{2\pi}} \exp \left(- \frac{x^{2}}{2} \right) dx \end{align}$

で求めることができます。 $\exp \left(- \frac{x^{2}}{2} \right)$ は偶関数であるため、

$\begin{align} \displaystyle \int_{-\sqrt{y}}^{\sqrt{y}} \frac{1}{\sqrt{2\pi}} \exp \left(- \frac{x^{2}}{2} \right) dx &= 2 \int_{0}^{\sqrt{y}} \frac{1}{\sqrt{2\pi}} \exp \left(- \frac{x^{2}}{2} \right) dx \\ &= 2 \left( \int_{-\infty}^{\sqrt{y}} \frac{1}{\sqrt{2\pi}} \exp \left(- \frac{x^{2}}{2} \right) dx - \int_{-\infty}^{0} \frac{1}{\sqrt{2\pi}} \exp \left(- \frac{x^{2}}{2} \right) dx \right) \end{align}$

ここで、ガウス積分の公式を用いると、

$\begin{align} \displaystyle \int_{-\infty}^{0} \frac{1}{\sqrt{2\pi}} \exp \left(- \frac{x^{2}}{2} \right) dx &= \int_{0}^{\infty} \frac{1}{\sqrt{2\pi}} \exp \left(- \frac{x^{2}}{2} \right) dx \\ &= \frac{1}{\sqrt{2\pi}} \cdot \sqrt{\frac{\pi}{2}} \\ &= \frac{1}{2} \end{align}$

となるため、

$\begin{align} \displaystyle P(X^{2} \le y) &= 2 \left( \int_{-\infty}^{\sqrt{y}} \frac{1}{\sqrt{2\pi}} \exp \left(- \frac{x^{2}}{2} \right) dx - \int_{-\infty}^{0} \frac{1}{\sqrt{2\pi}} \exp \left(- \frac{x^{2}}{2} \right) dx \right) \\ &= 2 \left( \phi(\sqrt{y}) - \frac{1}{2} \right) \end{align}$

が得られます。ここで $\phi (y)$ は、標準正規分布の累積分布関数です。

密度関数

密度関数 $g(y)$ は、累積分布関数を微分したものです。

$\displaystyle \frac{d}{dy} \phi(y) = \frac{1}{\sqrt{2 \pi}} \exp \left( - \frac{y}{2} \right)$

を用いると、

$\begin{align} \displaystyle g(y) &= \frac{dP}{dy} \\ &= \frac{d}{dy} 2 \left( \phi(\sqrt{y}) - \frac{1}{2} \right) \\ &= 2 \cdot \frac{1}{2\sqrt{y}} \cdot \frac{d}{dy} \phi (\sqrt{y}) \\ &= \frac{1}{\sqrt{2 \pi y}} \exp \left( - \frac{y}{2} \right) \quad (y \ge 0) \end{align}$

となります。

期待値

期待値 $E(Y)$ は、

$\begin{align} \displaystyle E(Y) &= \int_{-\infty}^{\infty} y \cdot g(y) dy \\ &= \int_{0}^{\infty} y \cdot \frac{1}{\sqrt{2 \pi y}} \exp \left( - \frac{y}{2} \right) dy \end{align}$

で求めますが、 $y = x^{2}$ と $dy = 2x \cdot dx$ であることを利用し、 $x$ が $0 \le x \le \infty$ の域値であること（累積分布関数の算出で偶関数の性質を使うときに、 $x$ が負の領域を考慮して累積分布関数を2倍していること）に注意すると、

$\begin{align} \displaystyle E(Y) &= \int_{0}^{\infty} y \cdot \frac{1}{\sqrt{2 \pi y}} \exp \left( - \frac{y}{2} \right) dy \\ &= \sqrt{\frac{2}{\pi}} \int_{0}^{\infty} x^{2} \exp \left( - \frac{x^{2}}{2} \right) \\ &= \sqrt{\frac{2}{\pi}} \cdot \frac{1}{4} \sqrt{\frac{\pi}{\left( \frac{1}{2} \right)^{3} }} \\ &= 1 \end{align}$

となります（ガウス積分の公式を利用）。

分散

$\begin{align} \displaystyle E(Y^{2}) &= \int_{-\infty}^{\infty} y^{2} \cdot g(y) dy \\ &= \int_{0}^{\infty} y^{2} \cdot \frac{1}{\sqrt{2 \pi y}} \exp \left( - \frac{y}{2} \right) dy \\ &= \sqrt{\frac{2}{\pi}} \int_{0}^{\infty} x^{4} \exp \left( - \frac{x^{2}}{2} \right) dx \\ &= \sqrt{\frac{2}{\pi}} \frac{3}{8} \sqrt{\frac{\pi}{\left( \frac{1}{2} \right)^{5} }} \\ &= 3 \end{align}$

となるため、分散 $V(Y)$ は、

$\begin{align} \displaystyle V(Y) &= E(Y^{2}) - E(Y)^{2} \\ &= 3 - 1 \\ &= 2 \end{align}$

と求まります。

5.8

コインの投げた回数を $x$ 回として、表が出た回数について $x \to a$ の極限を取ることを考えます。

右極限 $x \to a + 0$ を取ると $P(X = a)$ となります。ただし、累積分布関数は $G(a) = P(X \lt a)$ と定義されるため、 $\displaystyle \lim_{x \to a + 0} G(a) = P(X \le a) \neq P(X \lt a)$ より、右連続ではありません。

一方で左極限 $x \to a - 0$ を取ると $P(X = a - 1)$ となります。ここで、累積分布関数は $G(a) = P(X \lt a)$ と定義されるため、 $\displaystyle \lim_{x \to a - 0} G(a) = P(X \le a-1) = P(X \lt a)$ より、左連続になります（コインを投げた回数が離散値を取ることに注意します）。

一方、累積分布関数を $F(x) = P(X \le x)$ と定義した場合は、同様の議論から左連続ではなく右連続になります。

以上のことから、 $G(x)$ は右連続ではなく左連続になる点で $F(x)$ とは異なります。

2020-08-29

【統計学入門（東京大学出版会）】第4章練習問題解答

統計学

東京大学出版会から出版されている統計学入門（基礎統計学Ⅰ）について第4章の練習問題の解答を書いていきます。

本章以外の解答
4.1
- (i)
- (ii)
4.2
4.3
4.4
4.5
4.6
4.7
- (i)
- (ii)

本章以外の解答

本章以外の練習問題の解答は別の記事で公開しています。
必要に応じて参照してください。

4.1

(i)

6以外の目が出る確率は5/6であるため、さいころを4回投げたときに6の目が1回も出ない確率は、

$\displaystyle \left( \frac{5}{6} \right)^{4} = 0.482253$

となります。したがって、6の目が少なくとも1回出るほうに賭けたほうがよいことになります。

(ii)

さいころを2個同時に投げて(6, 6)の以外の目が出る確率は1/36であるため、24回投げたときに(6, 6)の目が1回も出ない確率は、

$\displaystyle \left( \frac{35}{36} \right)^{24} = 0.508596$

となります。したがって、(6, 6)の目が1回も出ないほうに賭けたほうがよいことになります。

4.2

さいころを2個同時に投げて目の和が12となる組合せは、(6, 6)の目しかありません。(6, 6)以外の目が出る確率は35/36です。
さいころを投げた回数を $n$ としたときに(6, 6)の目が出る確率が0.9を超えることは、次の式によって表すことができます。

$\displaystyle 1 - \left( \frac{35}{36} \right)^{n} > 0.9$

この式を満たす $n$ を求めると $n \gt 81.73637$ となるため、答えは82となります。なお、この不等式を求めるときに対数を取って計算すると楽ですが、符号の扱いに気を付けましょう。

4.3

30人から15人のグループを2つ作ることになるから、30人から15人を選ぶ組合せ数は次の式によって求められます。

$\displaystyle {}_{30}C_{15} = \frac{30!}{15! 15!}$

15と30の階乗は、スターリングの公式によって求めます。

$\displaystyle n! = \exp \left( n \log_{e} n - n + \frac{1}{2}\log_{e}2\pi n \right)$

30の階乗は $n=30$ として、次のように値を求めることができます。

$\displaystyle 30! = \exp \left( 30 \log_{e} 30 - 30 + \frac{1}{2}\log_{e}60\pi \right)$

同様に15の階乗は $n=15$ として、次のように値を求めます。

$\displaystyle 15! = \exp \left( 15 \log_{e} 15 - 15 + \frac{1}{2}\log_{e}30\pi \right)$

これらを計算すると、最終的に求めたい組合せ数 156415325.9563033 ≒ 156,415,326 が得られます。
ちなみに、Pythonプログラムを使うと簡単に階乗を求めることができます。

from math import factorial

combination = factorial(30) / factorial(15)**2
print(combination)

このプログラムを実行すると155,117,520が得られ、スターリングの公式を用いて求めた結果と実際の値とが近い値になっていることがわかります。

4.4

r人すべてが異なる誕生日である確率は、次の式で表されます。

$\displaystyle \prod_{i=1}^{r} \left( \frac{365 - (i -1)}{365} \right) = \prod_{i=1}^{r} \left( 1 - \frac{i - 1}{365} \right)$

r人のうち少なくとも2人が同じ誕生日である場合は、r人すべてが異なる誕生日である場合の補集合であるため、次の式で表されます。

$\displaystyle 1 - \prod_{i=1}^{r} \left( 1 - \frac{i - 1}{365} \right) \\ \displaystyle = 1 - 1 \times \left( 1 - \frac{1}{365} \right) \times \left( 1 - \frac{2}{365} \right) \times ... \times \left( 1 - \frac{r - 1}{365} \right) \\ \displaystyle = 1 - \left( 1 - \frac{1}{365} \right) \times \left( 1 - \frac{2}{365} \right) \times ... \times \left( 1 - \frac{r - 1}{365} \right)$

rが5, 10, 15, 20～25, 30, 35, 40, 50, 60のときの確率を求めるPythonプログラムを次に示します。

for rmax in [5, 10, 15, 20, 21, 22, 23, 24, 25, 30, 35, 40, 50, 60]:
    ans = 1
    for r in range(1, rmax+1):
        ans *= (1 - (r-1)/365)
    ans = 1 - ans
    print(f"r={rmax}, ans={ans}")

上記のプログラムを実行すると、次に示す結果が得られます。

r=5, ans=0.02713557369979347
r=10, ans=0.11694817771107768
r=15, ans=0.25290131976368646
r=20, ans=0.41143838358058016
r=21, ans=0.4436883351652059
r=22, ans=0.4756953076625502
r=23, ans=0.5072972343239854
r=24, ans=0.5383442579145288
r=25, ans=0.5686997039694639
r=30, ans=0.7063162427192687
r=35, ans=0.8143832388747152
r=40, ans=0.891231809817949
r=50, ans=0.9703735795779884
r=60, ans=0.994122660865348

4.5

ある人を起点にして考えると、他の人が隣に来る確率は2/3であるのに対して、向かい合いになる確率は1/3です。
問題の「向かい合って座るのが15組に対して、隣同士になるのが30組」になることは、この確率に沿っているだけと考えられます。したがって、心理学者の推論は妥当ではないといえます。

4.6

さいころの目の和が9になる場合と10になる場合のパターンについて、それぞれ発生しうるさいころの組合せ数を次に示します。

さいころの和が9

パターン	組合せ数
(1, 2, 6)	6
(1, 3, 5)	6
(1, 4, 4)	3
(2, 2, 5)	3
(2, 3, 4)	6
(3, 3, 3)	1

したがって、さいころの和が9になる組合せ数の合計は25となります。

さいころの和が10

パターン	組合せ数
(1, 3, 6)	6
(1, 4, 5)	6
(2, 2, 6)	3
(2, 3, 5)	6
(2, 4, 4)	3
(3, 3, 4)	3

したがって、さいころの和が10になる組合せ数の合計は27となります。

以上のことから、組合せ数が多いさいころの和が10のなる場合のほうが起こりやすいといえます。

4.7

(i)

ベイズの定理と全確率の法則から、

$\displaystyle P(C|A) = \frac{P(C) P(A|C)}{P(A)} = \frac{P(C) P(A|C)}{P(C) P(A|C) + P(C^{c}) P(A|C^{c})}$

$P(C^{c}) = 1 - P(C)$ 、 $P(A|C^{c}) = 1 - P(A^{c}|C^{c})$ であるから、

$\displaystyle P(C|A) = \frac{P(C) P(A|C)}{P(C) P(A|C) + P(C^{c}) (1 - P(A^{c}|C^{c})} = \frac{0.005 \times 0.95}{0.005 \times 0.95 + 0.995 \times (1 - 0.95)}$

となります。したがって、 $P(C|A) = 0.0872$ となります。

(ii)

$P(A|C) = P(A^{c}|C^{c}) = R$ として、

$\displaystyle P(C|A) = \frac{P(C) P(A|C)}{P(C) P(A|C) + P(C^{c}) (1 - P(A^{c}|C^{c})} = \frac{0.005 R}{0.005 R + 0.995 (1 - R)} \ge 0.90$

を満たす $R$ の範囲を求めると、 $R \ge 0.9994$ となります。

2020-08-22

【統計学入門（東京大学出版会）】第3章練習問題解答

統計学

前回に引き続き、東京大学出版会から出版されている統計学入門（基礎統計学Ⅰ）について第3章の練習問題の解答を書いていきます。

本章以外の解答
3.1
3.3
3.4
- (i)
- (ii)
- (iii)

本章以外の解答

本章以外の練習問題の解答は別の記事で公開しています。
必要に応じて参照してください。

3.1

自民得票率（tokuhyo）をX軸に、持ち家比率（motiie）をY軸にした散布図を描くPythonプログラムを次に示します。
matplotlibの図の背景はデフォルトで透明であるため、fig.patch.set_alpha(1) で不透明にしている点に注意してください。

import numpy as np
import matplotlib.pyplot as plt

tokuhyo = np.array([
    41.4, 76.3, 59.2, 51.8, 52.5, 53.2, 62.4, 55.0, 57.7,
    63.2, 37.5, 48.5, 32.4, 20.5, 47.9, 68.9, 68.5, 52.5, 63.3, 58.8,
    59.7, 48.4, 40.7, 51.0, 50.9, 34.3, 25.8, 32.1, 34.4, 55.1, 60.3,
    57.0, 45.6, 54.2, 55.1, 55.7, 70.3, 61.8, 47.6, 42.5, 71.3, 55.2,
    65.2, 42.9, 54.7, 62.0, 48.2
])
motiie = np.array([
    52.8, 71.2, 72.6, 63.7, 81.3, 81.8, 70.9, 74.0, 73.2,
    72.9, 66.7, 65.7, 43.7, 55.5, 79.6, 85.7, 75.3, 80.5, 73.0, 77.0,
    77.5, 69.2, 60.0, 78.2, 79.5, 61.8, 49.6, 59.6, 72.1, 71.0, 76.3,
    72.8, 71.8, 60.7, 67.0, 71.8, 71.2, 68.3, 68.5, 54.8, 76.0, 65.8,
    69.4, 66.9, 69.7, 71.2, 59.6
])

fig = plt.figure(figsize=(6, 6))
fig.patch.set_alpha(1)

plt.subplot(1, 1, 1)
plt.scatter(tokuhyo, motiie)
plt.xlabel("tokuhyo")
plt.ylabel("motiie")
plt.title("3.1")

上記のプログラムを実行すると、次の散布図が描画されます。

f:id:nuka137:20200822151604p:plain

変数 $x, y$ の間の相関係数 $r_{xy}$ は、次の式によって求められます。

$\displaystyle r_{xy} = \frac{\sum (x_{i} - \bar{x})(y_{i} - \bar{y})}{\sqrt{\sum (x_{i} - \bar{x})^{2}}\sqrt{\sum (y_{i} - \bar{y})^{2}}}$

変数 $x$ を自民得票率、変数 $y$ を持ち家比率として相関係数 $r_{xy}$ を求めるPythonプログラムを次に示します。

import numpy as np

tokuhyo = np.array([
    41.4, 76.3, 59.2, 51.8, 52.5, 53.2, 62.4, 55.0, 57.7,
    63.2, 37.5, 48.5, 32.4, 20.5, 47.9, 68.9, 68.5, 52.5, 63.3, 58.8,
    59.7, 48.4, 40.7, 51.0, 50.9, 34.3, 25.8, 32.1, 34.4, 55.1, 60.3,
    57.0, 45.6, 54.2, 55.1, 55.7, 70.3, 61.8, 47.6, 42.5, 71.3, 55.2,
    65.2, 42.9, 54.7, 62.0, 48.2
])
motiie = np.array([
    52.8, 71.2, 72.6, 63.7, 81.3, 81.8, 70.9, 74.0, 73.2,
    72.9, 66.7, 65.7, 43.7, 55.5, 79.6, 85.7, 75.3, 80.5, 73.0, 77.0,
    77.5, 69.2, 60.0, 78.2, 79.5, 61.8, 49.6, 59.6, 72.1, 71.0, 76.3,
    72.8, 71.8, 60.7, 67.0, 71.8, 71.2, 68.3, 68.5, 54.8, 76.0, 65.8,
    69.4, 66.9, 69.7, 71.2, 59.6
])

n = tokuhyo.size
tokuhyo_ave = np.mean(tokuhyo)
motiie_ave = np.mean(motiie)
rxy = np.sum((tokuhyo - tokuhyo_ave) * (motiie - motiie_ave)) / (np.sqrt(np.sum((tokuhyo - tokuhyo_ave)**2)) * np.sqrt(np.sum((motiie - motiie_ave)**2)))
print(f"Correlation Coefficient: {rxy}")

上記プログラムの実行結果を次に示します。

Correlation Coefficient: 0.6387815655787518

3.3

スピアマンの順位相関係数 $r_{s}$ は、次式で求められます。

$\displaystyle r_{s} = 1 - \frac{6}{n^{3} -n} \sum_{i = 1}^{n} (R_{i} - R_{i}^{'})^{2}$

一方、ケンドールの順位相関係数 $r_{k}$ は、次式で求められます。

$\displaystyle r_{k} = \frac{G - H}{n (n - 1) / 2}$

$G$ は、 $R_{i} \gt R_{j}$ かつ $R_{i}^{'} \gt R_{j}^{'}$ または、 $R_{i} \lt R_{j}$ かつ $R_{i}^{'} \lt R_{j}^{'}$ を満たす個数。
$H$ は、 $R_{i} \gt R_{j}$ かつ $R_{i}^{'} \lt R_{j}^{'}$ または、 $R_{i} \lt R_{j}$ かつ $R_{i}^{'} \gt R_{j}^{'}$ を満たす個数。

これらの式に従って $r_{s}$ と $r_{k}$ を求めるPythonプログラムを次に示します。

import numpy as np
import itertools


def spearman(d1, d2):
    n = d1.size
    return 1 - 6 * np.sum((d1 - d2)**2) / (n**3 - n)

def kendall(d1, d2):
    G = H = 0
    n = d1.size
    for i in range(n):
        Ri = d1[i]
        Rid = d2[i]
        for j in range(i+1, n):
            Rj = d1[j]
            Rjd = d2[j]
            if (Ri > Rj and Rid > Rjd) or (Ri < Rj and Rid < Rjd):
                G += 1
            elif (Ri > Rj and Rid < Rjd) or (Ri < Rj and Rid > Rjd):
                H += 1

    return (G - H) / (n*(n-1) / 2)

data = {
    "1": np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]),
    "2": np.array([1, 5, 2, 3, 6, 7, 15, 8, 4, 11, 10, 14, 18, 13, 22, 24, 16, 19, 30, 9, 25, 17, 26, 23, 12, 20, 28, 21, 27, 29]),
    "3": np.array([8, 3, 1, 4, 2, 5, 11, 7, 15, 9, 6, 13, 10, 22, 12, 14, 18, 19, 17, 22, 16, 24, 21, 20, 28, 30, 25, 26, 27, 29]),
    "4": np.array([20, 1, 4, 2, 6, 3, 12, 17, 8, 5, 18, 13, 23, 26, 29, 15, 16, 9, 10, 11, 30, 7, 27, 19, 14, 21, 28, 24, 22, 25]),
}

for d in itertools.combinations(data.keys(), 2):
    entry1 = d[0]
    entry2 = d[1]
    rs = spearman(data[entry1], data[entry2])
    rk = kendall(data[entry1], data[entry2])
    print(f"rs ({entry1} - {entry2}): {rs}")
    print(f"rk ({entry1} - {entry2}): {rk}")
    print("")

上記プログラムの実行結果を次に示します。

rs (1 - 2): 0.821579532814238
rk (1 - 2): 0.664367816091954

rs (1 - 3): 0.9272525027808676
rk (1 - 3): 0.767816091954023

rs (1 - 4): 0.5933259176863181
rk (1 - 4): 0.43908045977011495

rs (2 - 3): 0.671857619577308
rk (2 - 3): 0.5011494252873563

rs (2 - 4): 0.6373748609566184
rk (2 - 4): 0.46206896551724136

rs (3 - 4): 0.5343715239154616
rk (3 - 4): 0.3586206896551724

3.4

それぞれPythonのプログラムで解答を書きます。

(i)

乱数を発生させるPythonプログラムを示します。
1から11の乱数を発生させるために、random.randint を利用しました。

import random

random.seed(1)
d = int(random.randint(1, 11))

(ii)

変数 $x, y$ の間の相関係数 $r_{xy}$ は、次の式によって求められます。

$\displaystyle r_{xy} = \frac{\sum (x_{i} - \bar{x})(y_{i} - \bar{y})}{\sqrt{\sum (x_{i} - \bar{x})^{2}}\sqrt{\sum (y_{i} - \bar{y})^{2}}}$

相関係数 $r_{xy}$ を求めるPythonプログラムを次に示します。

import random
import numpy as np

random.seed(1)

x = np.array([71, 68, 66, 67, 70, 71, 70, 73, 72, 65, 66])
y = np.array([69, 64, 65, 63, 65, 62, 65, 64, 66, 59, 62])
x_mean = np.mean(x)
y_mean = np.mean(y)

x_sampled = np.zeros_like(x)
y_sampled = np.zeros_like(y)
for i in range(11):
    d = int(random.randint(1, 11)) - 1
    x_sampled[i] = x[d]
    y_sampled[i] = y[d]

rxy = np.sum((x_sampled - x_mean) * (y_sampled - y_mean)) / \
    (np.sqrt(np.sum((x_sampled - x_mean)**2)) * np.sqrt(np.sum((y_sampled - y_mean)**2)))
print(rxy)

上記プログラムの実行結果を次に示します。

0.49543369430686224

(iii)

(ii)を200回行い、 $r_{xy}$ の値をビン数10でヒストグラムにするPythonプログラムを次に示します。

import random
import numpy as np
import matplotlib.pyplot as plt

random.seed(1)

x = np.array([71, 68, 66, 67, 70, 71, 70, 73, 72, 65, 66])
y = np.array([69, 64, 65, 63, 65, 62, 65, 64, 66, 59, 62])
x_mean = np.mean(x)
y_mean = np.mean(y)

rxy_list = []
for _ in range(200):
    x_sampled = np.zeros_like(x)
    y_sampled = np.zeros_like(y)
    for i in range(11):
        d = int(random.randint(1, 11)) - 1
        x_sampled[i] = x[d]
        y_sampled[i] = y[d]

    rxy = np.sum((x_sampled - x_mean) * (y_sampled - y_mean)) / \
        (np.sqrt(np.sum((x_sampled - x_mean)**2)) * np.sqrt(np.sum((y_sampled - y_mean)**2)))
    rxy_list.append(rxy)

fig = plt.figure(figsize=(6, 6))
fig.patch.set_alpha(1)

plt.subplot(1, 1, 1)
plt.hist(rxy_list, bins=10)
plt.xlabel("rxy")
plt.ylabel("count")
plt.title("3.4 (iii)")

上記のプログラムを実行すると、次のような図が表示されます。

f:id:nuka137:20200822184620p:plain

2020-08-20

【統計学入門（東京大学出版会）】第2章練習問題解答

統計学

機械学習や深層学習を勉強する中で、その基礎となる統計学を抑えていないのはよくないと思い、東京大学出版会から出版されている統計学入門（基礎統計学Ⅰ）を読んでいます。執筆時点では第10章の練習問題を解いている途中ですが、本書は間違いなく良書だと思います。体系的に統計情報を学べるのはもちろんのこと、各章の最後にある練習問題はその章の内容を理解していないと解けない問題が多いため、手を動かすことでより深く理解できます。問題の中には、解いた結果の意味について考察するものもあり、本書の中で紹介された内容がどのように現実問題で役立つかを理解できます。自分は本書のおかげで、χの二乗分布が何に役立つのか理解できました。

しかしこの本、本の名前に「入門」と書いてはいますが、難易度は高めであることに注意が必要です。自分は大学が物理学科専攻で数式と格闘する毎日であったおかげか、数式で詰まることはほとんどありませんでしたが、大学教養レベルの数学の知識がないと読むのは厳しいと思います。

章末の練習問題の解答は、問題によっては答えのみしか書かれていないため不親切だと思いました。せっかく良問が揃っている中で解説がないのは、初学者にとって厳しいです。そこで、あとで問題を解き直したときに困らないように、練習問題の解答の詳解を記事としてまとめることにしました。第1章の練習問題など、自由研究のような問題を除いて解答に至るまでの途中の計算式などを書いていこうと思います。また、電卓を用いて値を求めるような問題に関しては、Pythonでプログラムを書いて求めています。

本章以外の解答
2.2
2.3
2.4

本章以外の解答

本章以外の練習問題の解答は別の記事で公開しています。
必要に応じて参照してください。

2.2

データA, B, Cを、平均差とジニ係数の式にそれぞれ代入すれば求まります。

平均差

$\displaystyle \sum_{i} \sum_{j} \frac{|x_{i} - x_{j}|}{n^{2}}$

平均差を求めるPythonプログラムを示します。

import numpy as np

A = np.array([0, 3, 3, 5, 5, 5, 5, 7, 7, 10])
B = np.array([0, 1, 2, 3, 5, 5, 7, 8, 9, 10])
C = np.array([3, 4, 4, 5, 5, 5, 5, 6, 6, 7])

def ave_diff(X):
    X_sum = 0
    for xi in X:
        for xj in X:
            X_sum += np.abs(xi - xj)
    return X_sum / (X.size**2)

print("AveDiff A: {}".format(ave_diff(A)))
print("AveDiff B: {}".format(ave_diff(B)))
print("AveDiff C: {}".format(ave_diff(C)))

上記プログラムの実行結果を次に示します。

AveDiff A: 2.76
AveDiff B: 3.76
AveDiff C: 1.2

ジニ係数

$\displaystyle GI = \sum_{i} \sum_{j} \frac{|x_{i} - x_{j}|}{2n^{2}\overline{x}}$

ジニ係数を求めるPythonプログラムを示します。

import numpy as np

A = np.array([0, 3, 3, 5, 5, 5, 5, 7, 7, 10])
B = np.array([0, 1, 2, 3, 5, 5, 7, 8, 9, 10])
C = np.array([3, 4, 4, 5, 5, 5, 5, 6, 6, 7])


def gini(X):
    X_mean = np.mean(X)
    X_sum = 0
    for xi in X:
        for xj in X:
            X_sum += np.abs(xi - xj)
    return X_sum / (2 * X.size**2 * X_mean)

print("Gini A: {}".format(gini(A)))
print("Gini B: {}".format(gini(B)))
print("Gini C: {}".format(gini(C)))

上記プログラムの実行結果を次に示します。

Gini A: 0.276
Gini B: 0.376
Gini C: 0.12

2.3

10年前と本年における学生の出身地の分布データを、エントロピー $H$ の式に代入して値を確認します。

$\displaystyle H(X) = - \sum_{i = 1}^{k} \hat{p_{i}} \log \hat{p_{i}}$

（ $p_{i}$ は、確率変数が $X = i$ となる確率です）

10年前と本年における学生の出身地の分布データについて、エントロピーを求めるPythonプログラムを示します。

import numpy as np

this_year = [32, 19, 10, 24, 15]
ten_years_ago = [28, 13, 18, 29, 12]


def entropy(X):
    f = np.array(X)
    p = f / np.sum(f)
    H = -np.sum(p * np.log10(p))
    return H

print("this year: {}".format(entropy(this_year)))
print("10 years ago: {}".format(entropy(ten_years_ago)))

上記プログラムの実行結果を次に示します。

this year: 0.667724435887455
10 years ago: 0.6704368955892825

エントロピーは確率変数（学生の出身地）の集中性を決めるため、10年前と本年では学生の出身地の集中性はそれほど変化していないと言えます。

2.4

データBを、標準得点と偏差値得点の式に代入します。

標準得点

$\displaystyle \mu = \sum^{n} x_{i} \\ \displaystyle \sigma = \sqrt{\sum^{n}(x_{i}^{2} - \mu^{2})} \\ \displaystyle z_{i} = \frac{x_{i} - \mu}{\sigma}$

$\mu$ は平均、 $\sigma$ は標準偏差、 $z_{i}$ は標準得点になります。

偏差値得点

$T_{i} = 10 z_{i} + 50$

$z_{i}$ は標準得点、 $T_{i}$ は偏差値得点になります。

これらの式を使って、データBの標準得点と偏差値得点を求めるプログラムを次に示します。

import numpy as np

B = np.array([0, 1, 2, 3, 5, 5, 7, 8, 9, 10])

B_ave = np.mean(B)
B_std = np.sqrt(np.sum(np.mean(B**2) - B_ave**2))
z = (B - B_ave) / B_std

print(f"Standard Score: {z}")
print("Deviation Score: {}".format(10 * z + 50))

上記プログラムの実行結果を次に示します。

Standard Score: [-1.52145155 -1.21716124 -0.91287093 -0.60858062  0.          0.
  0.60858062  0.91287093  1.21716124  1.52145155]
Deviation Score: [34.78548451 37.82838761 40.87129071 43.91419381 50.         50.
 56.08580619 59.12870929 62.17161239 65.21451549]

2020-08-18

Qiitaからはてなブログへ移行しました

その他

これまで技術関連の記事をQiitaに投稿していましたが、今後ははてなブログを使っていきたいと思います。 Qiitaに投稿した記事と同様、機械学習関連の技術記事を投稿していく予定です。

記念すべき最初の記事では、Qiitaの移行先としてはてなブログを選んだ理由について書きます。

Qiitaへの記事投稿をやめた理由
ブログサービスの比較
まとめ

Qiitaへの記事投稿をやめた理由

QiitaはSEOに対して強く、エンジニアで利用されている方も多いため、投稿した記事を多くのエンジニアから見てもらえる可能性が高いというメリットがありました。実際に私がQiitaに記事を投稿していたときも多くの方からLGTMをもらうことができ、そのたびに執筆のモチベーションが上がりました。

一方で最近では技術記事とは呼べない記事や、読まれることを前提としていないメモ書き程度の記事もよく見かけるようになりました（もちろん、素晴らしい技術記事を書かれている方は多数いますのであくまで一例です）。このような背景もあってQiitaに純粋な技術記事を投稿する意欲がなくなり、Qiitaへの記事投稿をやめて別のブログサービスへの移行を決めました。

ブログサービスの比較

noteやMediumなど個人でブログを開設する方法は多岐にわたるため、Qiitaの移行先としてどのブログサービスを選ぶか非常に悩みました。ここでは、著名なブログサービスの特徴を調べて理解した内容と、はてなブログを選ぶことになった理由を書きたいと思います。

note

手軽に記事を書けるところが良いと評判なブログサービスです。ユーザのつながりが強く、記事を見てもらいやすいところが利点であると、他の方が書かれたブログサービスの比較記事に書かれていました。 noteを本格的に利用していないため確かなことは言えませんが、noteを選ぶ1つのきっかけになるかもしれません。

また、noteでは有料記事や投げ銭システムなど、記事の収益化に便利な機能が提供されています。複数の記事をまとめることができるマガジン機能も有料化に対応しており、連載記事をまとめて販売するのもおもしろそうです。 noteの記事収益化に興味がある一方で、自分のような無名ユーザがいきなり有料記事を書いて収益を上げるのは難しそうだと感じました。

一方で技術系ブログを執筆するには、エディタが貧弱で厳しいものがありました。辛うじてソースコードのシンタックスハイライト機能があるものの、見出しレベルが1つしか設定できなかったり、リストが使えなかったり・・・といろいろなところで制約があります。機械学習関係の記事を書くとなると数式も多用することになるため、数式を画像化してアップロードして・・・という作業が繰り返されることになると思うと地獄を見られそうです。

Medium

海外のユーザが多いブログサービスです。海外ユーザが多いことから、英語で記事を執筆することでより多くの人に記事を見てもらえるというコメントも見られました。

海外の企業が公式のブログとしてMediumを利用していることも多く、品質の高い記事を読むことができる点も特徴でしょう。著名な機械学習ライブラリPyTorchもMediumを使っていますし、独自のブログに移行する前はTensorFlowもMediumを利用していました。

Mediumは、clap（QiitaのLGTMに相当）の数に対して収益が得られる仕組みが面白いと思いました。 QiitaでLGTMをもらえるだけでも記事を執筆するモチベーションが上がりましたが、それに加えて収益も得られるとなるとさらにモチベーションが上がりそうです。

WordPress

WordPressは、ブログの枠組みを超えてさまざまなコンテンツを管理できるソフトウェアです。

カスタマイズの自由度が高いことがWordPressのメリットです。一方で、WordPressを動作させるためのレンタルサーバやWordPress自体のメンテナンスに時間を割かれることがデメリットになります。 WordPressをメンテナンスするコストのほうが、記事を書くコストよりも大きくなる可能性があるため、ある程度本格的にWebコンテンツを作ることを想定していない人にとっては導入の敷居が高いように思えます。

はてなブログ

はてなブログは、はてなダイアリーが前身ということもあり、技術系のブログとして採用されることが多い印象です。

技術系のブログを採用するにあたってエディタは大切な要素ですが、技術系のブログが多いはてなブログではあまりその心配はなさそうです。また、（個人的な印象で）はてなブログは収益化が難しいと感じていましたが、最近広告による収入も可能であることを知り、ブログのアクセス数が増えれば収益化も不可能ではないと感じました。

まとめ

代表的なブログサービスについて、その特徴をまとめると次のようになります。

ブログサービス	特徴
note	・記事投稿までの流れがシンプル・記事収益化の手続きが簡単・ユーザ間のつながりが強い（？）・エディタが貧弱
Medium	・海外ユーザが多い・企業の公式ブログとして利用されている・clapによる収益化が可能・日本では知名度が低い
WordPress	・カスタマイズの幅が広い・メンテナンスコストが高い
はてなブログ	・技術系のブログ多い・広告による収益化が可能

自分は次の理由から「はてなブログ」を選びました。

技術系の記事を中心に投稿する予定である（技術系ブログとしての実績やエディタが重要）
ブログ開設時点では個人としての知名度が低く、有料記事として利益を上げられる見込みがない
日本国内を対象としている
現時点でブログの収益化は考えていないが、知名度が高くなったときのことも踏まえて収益化できる仕組みが欲しい

なお、Qiitaに投稿した記事はそのまま残しますが、今後ははてなブログに技術記事を投稿していこうと思います。もし投稿した記事が役立った方は、スターやブクマしてもらえると執筆のモチベーションが上がります！