尤度・最尤推定

2020-03-01 / tau / コメントする

コインの例

最尤推定法(maximum likelihood estimation)は確率分布 $f_D$ のパラメーター $\theta$ を、その母集団から得られたサンプルから推定する方法。

たとえば、表の出る確率が必ずしも1/2でないコインを50回投げた時、表が15回出たとすると、このコインの表が出る確率はいくらか、といった問題。20～30回くらいの間だと普通のコインで表／裏の確率は1/2かなと思うが、50回投げて表が15回しかでないとちょっと疑いたくなる。

ここで、表が出る確率を $p$ として、 $n$ 回の試行中表が $k$ 回出る確率は以下の通り。

(1) $\begin{equation*} P(k:n\,|\,p) = \binom{n}{k} p^k (1-p)^{n-k} \end{equation*}$

先の問題の結果について、母集団の確率をいくつか仮定して計算してみると

(2) $\begin{equation*} \begin{align} p &= 1/2 \quad \rightarrow \quad P(15:50\,|\,0.5) \approx 0.002 \\ p &= 1/3 \quad \rightarrow \quad P(15:50\,|\,0.3) \approx 0.107 \\ p &= 1/4 \quad \rightarrow \quad P(15:50\,|\,0.1) \approx 0.089 \end{align} \end{equation*}$

得られたサンプルの下では、少なくともぼ表が出る確率を1/2や1/4と考えるよりは、1/3と考えた方がこのサンプルのパターンが発生する確率が高い。

そこで、サンプルのパターンに対して母集団の確率を変化させてみると、以下のようになる。

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

n = 50
k = 15
p = np.linspace(0, 1, 100)

fig, ax = plt.subplots(figsize=(4.8, 3.6))

ax.plot(p, binom.pmf(k, n, p), c='navy')
ax.plot((k/n, k/n), (0, binom.pmf(k, n, k/n)), c='gray', linestyle='dashed')

ax.set_xlim(0, 1)
ax.set_ylim(0, 0.15)
ax.set_xticks(np.arange(0, 1.1, 0.1))

plt.show()

import numpy as np

import matplotlib.pyplot as plt

from scipy.stats import binom

n = 50

k = 15

p = np.linspace(0, 1, 100)

fig, ax = plt.subplots(figsize=(4.8, 3.6))

ax.plot(p, binom.pmf(k, n, p), c='navy')

ax.plot((k/n, k/n), (0, binom.pmf(k, n, k/n)), c='gray', linestyle='dashed')

ax.set_xlim(0, 1)

ax.set_ylim(0, 0.15)

ax.set_xticks(np.arange(0, 1.1, 0.1))

plt.show()

グラフ上は、母集団確率が15/50=0.3のときにサンプルパターンの発生確率が最も高くなっている。

解析的に確率が最大となる母集団確率を求めるために、式(1)を $p$ で微分してゼロとなる点を求める。

(3) $\begin{gather*} \frac{d P(k:n\,|\,p)}{dp} = \binom{n}{k}\left(kp^{k-1} (1-p)^{n-k} - p^k (n-k)(1-p)^{n-k-1} \right) = 0 \\ p^{k-1} (1-p)^{n-k-1} \left( k (1-p) - (n-k) p \right) = 0 \\ p^{k-1} (1-p)^{n-k-1} \left( k - np \right) = 0 \end{gather*}$

上式の解は $p = 0, 1, k/n$ となるが、 $p=0, 1$ の場合は式(1)がゼロとなることから、確率が最大となるのは $p=k/n$ であることがわかる。

尤度関数と最尤法

母集団の確率分布を $f$ 、母集団のパラメータを $\theta$ 、得られたサンプルセットを $x_1, \ldots, x_n$ とすると、パラメーターが $\theta$ である確率は

(4) $\begin{equation*} Pr(x_1, \ldots, x_n \,|\, \theta) = f(x_1, \ldots, x_n \,|\, \theta) \end{equation*}$

このとき、パラメーター $\theta$ に関する尤度関数は以下のように定義される。

(5) $\begin{equation*} L(\theta) = f(x_1, \ldots, x_n \,|\, \theta) \end{equation*}$

通常、尤度関数が最大となるパラメーターを求めるためには対数尤度関数が用いられる。

(6) $\begin{equation*} \frac{d}{d \theta} \ln L(\theta) = 0 \end{equation*}$

サンプルの質と尤度

先のコインの問題で、試行数が10回の場合と100回の場合を比べてみる。それぞれについて、表の確率が0.3、0.5、0.7に対応する回数は試行回数が10回の場合は3回、5回、7回、試行回数100回に対しては30回、50回、70回である。それぞれの表の階数に対して、母集団確率の尤度関数を描くと以下のようになる。

import numpy as np
import matplotlib.pyplot as plt
from functools import reduce

def comb(n, k):
    numer = reduce(lambda x,y:x*y, range(n - k + 1, n))
    denom = reduce(lambda x,y:x*y, range(1, k + 1))
    return numer // denom

def binom(n, k, p):
    return comb(n, k) * p**k * (1 - p)**(n - k)

n = [10, 100]
k = [np.array([3, 5, 7]), np.array([30, 50, 70])]

p = np.linspace(0, 1, 100)

fig, axs = plt.subplots(1, 2, figsize=(9.6, 3.6))

for ax, n, k in zip(axs, n, k):
    ax.set_title("n={}".format(n))

    for t in k:
        p_mle = t / n
        prob = binom(n, t, p)

        ax.plot(p, prob, c='navy')
        ax.plot((p_mle, p_mle), (0, binom(n, t, p_mle)), c='gray', linestyle='dashed')

        ax.set_xlim(0, 1)
        ax.set_xticks(np.arange(0, 1.1, 0.1))

plt.show()

import numpy as np

import matplotlib.pyplot as plt

from functools import reduce

def comb(n, k):

numer = reduce(lambda x,y:x*y, range(n - k + 1, n))

denom = reduce(lambda x,y:x*y, range(1, k + 1))

return numer // denom

def binom(n, k, p):

return comb(n, k) * p**k * (1 - p)**(n - k)

n = [10, 100]

k = [np.array([3, 5, 7]), np.array([30, 50, 70])]

p = np.linspace(0, 1, 100)

fig, axs = plt.subplots(1, 2, figsize=(9.6, 3.6))

for ax, n, k in zip(axs, n, k):

ax.set_title("n={}".format(n))

for t in k:

p_mle = t / n

prob = binom(n, t, p)

ax.plot(p, prob, c='navy')

ax.plot((p_mle, p_mle), (0, binom(n, t, p_mle)), c='gray', linestyle='dashed')

ax.set_xlim(0, 1)

ax.set_xticks(np.arange(0, 1.1, 0.1))

plt.show()

試行回数が100回の場合には分布が尖っており、最尤推定されたパラメーターの確度が高いが、10回の場合には分布がなだらかで、最尤推定された値以外のパラメーターの確率も高くなっている。

最尤推定の場合、サンプルの規模が推定値の信頼度に関わってくる。

二項分布の平均と分散

2019-12-07 / tau / コメントする

概要

二項分布 $B(n, p)$ の平均と分散は以下のようになる。

(1) $\begin{alignat*}{1} E(X) &= np \\ V(X) &= np(1-p) \end{alignat*}$

これらを導くのに、有用なテクニックを使っているのでまとめておく。

直接定義式から導く方法

この方法は、平均、分散の定義式から直接導いていく過程で、意図的に二項展開の形に持ち込んでいく方法。

平均

平均の定義から、以下のように変形していく。

(2) $\begin{alignat*}{1} E(X) &= \sum_{k=0}^n k {}_n\mathrm{C}_k p^k (1-p)^{n-k} \\ &= \sum_{k=1}^n k \frac{n!}{(n-k)! k!} p^k (1-p)^{n-k} \\ &= \sum_{k=1}^n n \frac{(n-1)!}{\left( (n-1) - (k-1) \right)! (k-1)!} pp^{k-1} (1-p)^{(n-1)-(k-1)} \end{alignat*}$

上式では、各項に $k$ が乗じられていることから、以下の流れで変形している。

$k=0$ のとき1番目の項はゼロとなるので、和の開始値を $k=1$ とする
分子にある $k$ を使って、 ${}_n C_k$ の分母において $k! \rightarrow (k-1)!$ とする
これに整合させるため、分母において $(n-k)!$ を $((n-1) - (k-1))!$ と変形
さらに組み合わせの式に整合させるため分子を $n (n-1)!$ と変形し、最終的に ${}_{n-1} C_{k-1}$ を引き出している。
後に二項定理を使うため、 $p, \; 1-p$ の指数も調整している

ここで $k-1 = k'$ とおくと、カウンターの範囲は $k=1 \sim n$ から $k'=0 \sim n-1$ となることから、

(3) $\begin{alignat*}{1} E(X) &= np \sum_{k'=0}^{n-1} \frac{(n-1)!}{\left( (n-1) - k' \right)! k'!} p^{k'} (1-p)^{(n-1)-k'} \\ &= np (p + 1 - p)^{n-1} \\ &= np \end{alignat*}$

上式では、変形した和の部分が二項展開の形になっていることを利用している。

分散

分散については、 $k^2$ が各項に乗じられるが、これを $k(k-1) + k$ と変形して、階乗のランクを下げるところがミソ。

(4) $\begin{alignat*}{1} V(X) &= E(X^2) - [E(X)]^2 \\ &= \sum_{k=0}^n k^2 {}_n C_k p^k (1-p)^{n-k} - (np)^2 \\ &= \sum_{k=0}^n \left( k(k-1) + k \right) {}_n C_k p^k (1-p)^{n-k} - (np)^2 \\ &= \sum_{k=2}^n k(k-1) {}_n C_k p^k (1-p)^{n-k} + \sum_{k=1}^n k {}_n C_k p^k (1-p)^{n-k} - (np)^2 \\ \end{alignat*}$

上式で、1項目は $k(k-1)$ が乗じられているのでカウンターを $k=2$ から、2項目は同じく $k=1$ からとしている。

ここで1項目について $k'' = k - 2$ と置いて、平均の時と同じ考え方で以下のように変形。

(5) $\begin{alignat*}{1} &\sum_{k=2}^n k(k-1) {}_n C_k p^k (1-p)^{n-k} \\ &= \sum_{k=2}^n k(k-1) \frac{n!}{(n-k)! k!} p^k (1-p)^{n-k} \\ &= \sum_{k=2}^n n(n-1) \frac{(n-2)!}{((n-2) -(k-2))! (k-2)!} p^2 p^{k-2} (1-p)^{(n-2)-(k-2)} \\ &= n(n-1) p^2 \sum_{k'=0}^{n-2} \frac{(n-2)!}{((n-2) -k'')! k''!} p^{k''} (1-p)^{(n-2)-k''} \\ &= n(n-1) p^2 ( p + (1-p))^{n-2} \\ &= n(n-1) p^2 \end{alignat*}$

2項目については、 $k' = k-1$ とおいて、

(6) $\begin{alignat*}{1} &\sum_{k=1}^n k {}_n C_k p^k (1-p)^{n-k} \\ &= \sum_{k=1}^n n \frac{(n-1)!}{\left( (n-1) - (k-1) \right) ! (k-1)!} p p^{k-1} (1-p)^{(n-1)-(k-1)} \\ &= np \sum_{k'=0}^{n-1} \frac{(n-1)!}{\left( (n-1) - k'\right) ! k'!} p p^{k'} (1-p)^{(n-1)-k'} \\ &= np \left( (p + (1-p) \right) ^{n-1} \\ &= np \end{alignat*}$

以上を併せて、

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

微分による方法

この方法は、 $kp^k, \; k^2 p^k$ の形に着目して、全事象の式を微分する方法。式展開が素直であり、平均の式を微分した結果がそのまま分散の式になってしまうところが美しい。

平均

二項分布の全確率の和は1となる。

(8) $\begin{equation*} \sum_{k=0}^n {}_n \mathrm{C}_k p^k (1-p)^{n-k} = (p + 1 - p)^n = 1 \end{equation*}$

この式の両辺を $p$ で微分する。

(9) $\begin{gather*} \sum_{k=0}^n {}_n \mathrm{C}_k \left( k p^{k-1} (1-p)^{n-k} - (n-k)p^k (1-p)^{n-k-1} \right) = 0 \\ \sum_{k=0}^n {}_n \mathrm{C}_k p^{k-1} (1-p)^{n-k-1} \left( k(1-p) - (n-k)p \right) = 0 \\ \sum_{k=0}^n {}_n \mathrm{C}_k p^{k-1} (1-p)^{n-k-1} (k - np) = 0 \\ \end{gather*}$

両辺に $p(1-p)$ をかける。

(10) $\begin{gather*} \sum_{k=0}^n {}_n \mathrm{C}_k p^k (1-p)^{n-k} (k - np) = 0 \\ \sum_{k=0}^n k {}_n \mathrm{C}_k p^k (1-p)^{n-k} = np \sum_{k=0}^n {}_n \mathrm{C}_k p^k (1-p)^{n-k}\\ \therefore E(X) = np \end{gather*}$

分散

式(10)をもう一度 $p$ で微分する。

(11) $\begin{gather*} \sum_{k=0}^n k {}_n \mathrm{C}_k p^k (1-p)^{n-k} = np \\ \sum_{k=0}^n k {}_n \mathrm{C}_k \left( kp^{k-1} (1-p)^{n-k} -(n-k)p^k (1-p)^{n-k-1} \right) = n\\ \sum_{k=0}^n k {}_n \mathrm{C}_k p^{k-1} (1-p)^{n-k-1} \left( k (1-p) -(n-k)p \right) = n\\ \sum_{k=0}^n k {}_n \mathrm{C}_k p^{k-1} (1-p)^{n-k-1} (k - np) = n\\ \end{gather*}$

両辺に $p(1-p)$ をかける。

(12) $\begin{gather*} \sum_{k=0}^n k {}_n \mathrm{C}_k p^k (1-p)^{n-k} (k - np) = np(1-p) \\ \sum_{k=0}^n k^2 {}_n \mathrm{C}_k p^k (1-p)^{n-k} - np \sum_{k=0}^n k {}_n \mathrm{C}_k p^k (1-p)^{n-k}= np(1-p) \\ E(X^2) - (np)^2 = np(1-p) \\ \therefore V(X) = np(1-p) \end{gather*}$

回帰分析～重回帰

2019-11-18 / tau / コメントする

概要

単回帰が $(X, Y)$ への線形関係を扱うのに対して、重回帰は複数の説明変数 $X_k \; (k=1 \ldots m)$ と $Y$ の線形関係を扱う。

データ $(X_1, \ldots , X_m, Y) = (x_{1i}, \ldots, x_{mi}, y_i), (i = 1 \ldots n)$ に対して以下の線形式で最も説明性の高いものを求める。

(1) $\begin{equation*} \hat{y} = w_0 + w_1 x_1 + \cdots + w_m x_m \end{equation*}$

そのために、データとその推測値の残差 $\hat{y_i} - y_i$ の平方和が最小となるような係数 $w_0, w_1, \ldots, w_m$ を最小二乗法により求める。

定式化

説明変数とターゲットのデータセット $n$ 組が次のように得られているとする。

(2) $\begin{equation*} x_{1i}, \ldots , x_{mi}, y_i \quad {\rm where} \quad i = 1, \ldots , n \end{equation*}$

残差の平方和については

(3) $\begin{equation*} \min \sum_{i=1}^n (\hat{y} -y_i)^2 \quad {\rm where} \quad \hat{y}_i = w_0 + w_1 x_{1i} + \cdots w_m x_{mi} \end{equation*}$

ここで残差を最小化する $w_0, w_1, \ldots, w_m$ を求めるために、それぞれで偏微分する。

(4) $\begin{equation*} \begin{array}{c} \left\{ \begin{array}{l} \displaystyle \frac{\partial}{\partial w_0} \sum_{i=1}^n (w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i)^2 = 0 \\ \displaystyle \frac{\partial}{\partial w_1} \sum_{i=1}^n (w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i)^2 = 0 \\ \vdots \\ \displaystyle \frac{\partial}{\partial w_m} \sum_{i=1}^n (w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i)^2 = 0 \\ \end{array} \right. \end{equation*}$

変形すると、

(5) $\begin{equation*} \left\{ \begin{array}{l} \displaystyle \sum_{i=1}^n 2(w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i) = 0 \\ \displaystyle \sum_{i=1}^n 2(w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i) x_{1i} = 0 \\ \vdots \\ \displaystyle \sum_{i=1}^n 2(w_0 + w_1 x_{1i} + \cdots + w_m x_{mi} - y_i) x_{mi} = 0 \\\end{array} \right. \end{array} \end{equation*}$

更にこれを展開するのに、 $S_k = \sum_{i=1}^n x_{ki}$ , $S_{kl} = \sum_{i=1}^n x_{ki} x_{li}$ , $S_y = \sum_{i=1}^n y_i$ , $S_{ky} = \sum_{i=1}^n x_{ki} y_i$ とおいて

(6) $\begin{equation*} \left\{ \begin{array}{l} n w_0 + w_1 S_1+ \cdots + w_m S_m - S_y = 0 \\ w_0S_1 + w_1 S_{11} + \cdots + w_m S_{1m} - S_{1y} = 0 \\ \vdots \\ w_0 S_m + w_1 S_{m1} + \cdots + w_m S_{mm} - S_{my} = 0 \\ \end{array} \right. \end{array} \end{equation*}$

これを行列形式で表示すると

(7) $\begin{equation*} \left( \begin{array}{cccc} n & S_1 & \cdots & S_m \\ S_1 & S_{11} & \cdots & S_{1m} \\ \vdots & \vdots & & \vdots \\ S_m & S_{m1} & \cdots & S_{mm} \end{array} \right) \left( \begin{array}{c} w_0 \\ w_1 \\ \vdots \\ w_m \end{array} \right) = \left( \begin{array}{c} S_y \\ S_{1y} \\ \vdots \\ S_{my} \end{array} \right) \end{equation*}$

解の導出

式(7)の連立方程式を解けば各係数を得ることができる。ここで左辺の係数行列は以下のようにも表せる。

(8) $\begin{equation*} \left( \begin{array}{cccc} 1 & 1 & \cdots & 1 \\ x_{11} & x_{12} & \cdots & x_{1n} \\ \vdots & \vdots & & \vdots \\ x_{m1} & x_{m2} & \cdots & x_{mn} \end{array} \right) \left( \begin{array}{cccc} 1 & x_{11} & \cdots & x_{m1} \\ 1 & x_{12} & \cdots & x_{m2} \\ \vdots & \vdots & & \vdots \\ 1 & x_{1n} & \cdots & x_{mn} \end{array} \right) = \boldsymbol{X}^T \boldsymbol{X} \end{equation*}$

また右辺は次のように表せる。

(9) $\begin{equation*} \left( \begin{array}{cccc} 1 & 1 & \cdots & 1 \\ x_{11} & x_{12} & \cdots & x_{1n} \\ \vdots & \vdots & & \vdots \\ x_{m1} & x_{m2} & \cdots & x_{mn} \end{array} \right) \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array} \right) = \boldsymbol{X}^T \boldsymbol{y} \end{equation*}$

したがって式(7)は以下のように表せる。

(10) $\begin{equation*} \boldsymbol{X}^T \boldsymbol{X} \boldsymbol{w} = \boldsymbol{X}^T \boldsymbol{y} \end{equation*}$

これを解くと

(11) $\begin{equation*} \boldsymbol{w} = ( \boldsymbol{X}^T \boldsymbol{X} )^{-1} \boldsymbol{X}^T \boldsymbol{y} \end{equation*}$

ベクトル・行列による表現と解

定式化の段階からベクトル・行列を用いてみる。式(1)の各変数を全てベクトル形式で表すと以下のようになる。

(12) $\begin{equation*} \hat{\boldsymbol{y}} = \boldsymbol{x}^T \boldsymbol{w} \end{equation*}$

ここで、

(13) $\begin{equation*} \boldsymbol{x} = \left( \begin{array}{c} 1 & x_1 & \vdots & x_n \end{array} \right) \quad , \quad \boldsymbol{w} = \left( \begin{array}{c} w_0 & w_1 & \vdots & w_n \end{array} \right) \end{equation*}$

複数データからなるデータセットに対して、式(8)の行列表現を用いると、以下のように表現できる。

(14) $\begin{equation*} \hat{\boldsymbol{y}} = \boldsymbol{X} \boldsymbol{w} \end{equation*}$

各データの残差については以下のようになる。

(15) $\begin{equation*} \hat{\boldsymbol{y}} - \boldsymbol{y} = \boldsymbol{X} \boldsymbol{w} - \boldsymbol{y} \end{equation*}$

ここで残差の二乗和を損失関数Lとして、これを行列で表現すると以下のようになる。

(16) $\begin{align*} L &= \left( \hat{\boldsymbol{y}} - \boldsymbol{y} \right)^T \left( \hat{\boldsymbol{y}} - \boldsymbol{y} \right) \\ &= \left( \boldsymbol{Xw} - \boldsymbol{y} \right)^T \left( \boldsymbol{Xw} - \boldsymbol{y} \right) \\ &= \left( \boldsymbol{w}^T \boldsymbol{X}^T - \boldsymbol{y}^T \right) \left( \boldsymbol{Xw} - \boldsymbol{y} \right) \\ &= \boldsymbol{w}^T \boldsymbol{X}^T \boldsymbol{Xw} - \boldsymbol{w}^T \boldsymbol{X}^T \boldsymbol{y} - \boldsymbol{y}^T \boldsymbol{Xw} + \boldsymbol{y}^T \boldsymbol{y} \end{align*}$

上式の第2項、第3項については以下のように整理できる。

(17) $\begin{align*} - \boldsymbol{w}^T \boldsymbol{X}^T \boldsymbol{y} - \boldsymbol{y}^T \boldsymbol{Xw} = - \left( \boldsymbol{Xw} \right)^T \boldsymbol{y} - \boldsymbol{y}^T \boldsymbol{Xw} = -2 \boldsymbol{y}^T \boldsymbol{Xw} \end{align*}$

損失関数Lはスカラーであり、これを最小にするためにベクトルwで微分し、それらの値がゼロとなるような方程式とする。

(18) $\begin{align*} \frac{dL}{d\boldsymbol{w}} = \frac{d}{d\boldsymbol{w}} \left( \boldsymbol{w}^T \boldsymbol{X}^T \boldsymbol{Xw} - 2 \boldsymbol{y}^T \boldsymbol{Xw} + \boldsymbol{y}^T \boldsymbol{y} \right) = 2 \boldsymbol{X}^T \boldsymbol{Xw} -2 \boldsymbol{X}^T \boldsymbol{y} = \boldsymbol{0} \end{align*}$

なお第1項については、式(7)と式(8)より $\boldsymbol{X}^T\boldsymbol{X}$ は対象行列なので転置しても同じ行列となり、行列の微分の公式から以下のようになることを利用している。

(19) $\begin{align*} \frac{d}{d\boldsymbol{x}} \boldsymbol{w}^T \boldsymbol{X}^T \boldsymbol{Xw} = \left( \boldsymbol{X}^T \boldsymbol{X} + \left(\boldsymbol{X}^T \boldsymbol{X} \right)^T \right) \boldsymbol{w} = 2 \boldsymbol{X}^T \boldsymbol{X} \boldsymbol{w} \end{align*}$

式(18)をwについて解いて以下を得る。

(20) $\begin{equation*} \boldsymbol{w} = \left( \boldsymbol{X}^T \boldsymbol{X} \right)^{-1} \boldsymbol{X}^T \boldsymbol{y} \end{equation*}$

多重共線性

式(7)から $w_0$ を消去する。

(21) $\begin{equation*} \left( \begin{array}{ccc} S_{11} - \dfrac{{S_1}^2}{n} & \cdots & S_{1m} - \dfrac{S_1 S_m}{n}\\ \vdots & & \vdots \\ S_{m1} - \dfrac{S_m S_1}{n} & \cdots & S_{mm} - \dfrac{{S_m_}^2}{n} \end{array} \right) \left( \begin{array}{c} w_1 \\ \vdots \\ w_m \end{array} \right) = \left( \begin{array}{c} S_{1y} - \dfrac{S_1 S_y}{n} \\ \vdots \\ S_{my} - \dfrac{S_m S_y}{n} \end{array} \right) \end{equation*}$

これを分散・共分散で表すと以下の通り。ただし以下で $V_{ii}=\rm{Var}(X_i)$ 、 $V_{ij} = \rm{Cov}(X_i, X_j)$ と表している。

(22) $\begin{equation*} \left( \begin{array}{ccc} V_{11} & \cdots & V_{1m}\\ \vdots & & \vdots \\ V_{m1} & \cdots & V_{mm} \end{array} \right) \left( \begin{array}{c} w_1 \\ \vdots \\ w_m \end{array} \right) = \left( \begin{array}{c} V_{1y} \\ \vdots \\ V_{my} \end{array} \right) \end{equation*}$

ここで $x_j = a x_i + b$ と2つの間に完全な線形関係がある場合、分散・共分散の性質から以下の関係が成り立つ。

(23) $\begin{equation*} V_{jj} = a^2V_{ii}, \; V_{ji} = V_{ij} = aV_{ii}, \; V_{jk} = V_{kj} = aV_{ji} = aV_{ij} \end{equation*}$

これらを式(22)の係数行列に適用すると、

(24) $\begin{align*} &\left[ \begin{array}{ccccccc} V_{11} & \cdots & V_{1i} & \cdots & V_{1j} & \cdots & V_{1m}\\ \vdots && \vdots && \vdots && \vdots\\ V_{i1} & \cdots & V_{ii} & \cdots & V_{ij} & \cdots & V_{im}\\ \vdots && \vdots && \vdots && \vdots\\ V_{j1} & \cdots & V_{ji} & \cdots & V_{jj} & \cdots & V_{jm}\\ \vdots && \vdots && \vdots && \vdots\\ V_{m1} & \cdots & V_{mi} & \cdots & V_{mj} & \cdots & V_{mm} \end{array} \right]\\ &= \left[ \begin{array}{ccccccc} V_{11} & \cdots & V_{1i} & \cdots & aV_{1i} & \cdots & V_{1m}\\ \vdots && \vdots && \vdots && \vdots\\ V_{i1} & \cdots & V_{ii} & \cdots & aV_{ii} & \cdots & V_{im}\\ \vdots && \vdots && \vdots && \vdots\\ aV_{i1} & \cdots & aV_{ii} & \cdots & a^2V_{ii} & \cdots & aV_{im}\\ \vdots && \vdots && \vdots && \vdots\\ V_{m1} & \cdots & V_{mi} & \cdots & aV_{mi} & \cdots & V_{mm} \end{array} \right] \end{align*}$

係数行列の行(列)が一次従属となっていて、その行列式がゼロとなることから、連立方程式は解を持たない。変数間の関係が線形でなくても相関がかなり高い場合は行列式がゼロに近づき、解が不安定になる。

単回帰分析～コンビニエンスストア

2019-09-29 / tau / コメントする

概要

Pythonを使った統計計算と図示の練習のため、コンビニエンスストアで単回帰分析をやってみた。

コンビニの店舗数は「商業動態統計年報」の2016年データを使い、説明変数として2015年の国勢調査人口、2018年の国土地理院による都道府県面積、2017年道路統計年報の道路実延長データを使った。

計算コードは以下の通り。

import numpy as np
import matplotlib.pyplot as plt

master_data = np.loadtxt('conv-store-data.csv', delimiter=',',
        skiprows=2, usecols=(1, 2), encoding='utf8')
print(master_data)

y = num_stores = master_data[:,0]
x = population = master_data[:,1]

sxy = np.sum((x - x.mean())*(y - y.mean()))
sxx = np.sum((x - x.mean())**2)
syy = np.sum((y - y.mean())**2)
a = sxy / sxx
b = y.mean() - a*x.mean()
r = sxy / np.sqrt(sxx * syy)
n = len(num_stores)
ESS = n * sxy**2 / sxx
TSS = n * syy
R2 = ESS / TSS

max_population = 15000000
max_stores = 8000
plt.xlim(0, max_population)
plt.ylim(0, max_stores)
plt.xlabel('population')
plt.ylabel('number of stores')

plt.scatter(population, num_stores, label='#stores')

x = np.linspace(0, max_population)
plt.plot(x, a*x + b, color='m')

plt.text(1000000, 7600, 'a = ' + str(a))
plt.text(1000000, 7300, 'b = ' + str(b))
plt.text(1000000, 7000, 'r = '+str(r))
plt.text(1000000, 6700, 'R2 = '+str(R2))

plt.show()

import numpy as np

import matplotlib.pyplot as plt

master_data = np.loadtxt('conv-store-data.csv', delimiter=',',

skiprows=2, usecols=(1, 2), encoding='utf8')

print(master_data)

y = num_stores = master_data[:,0]

x = population = master_data[:,1]

sxy = np.sum((x - x.mean())*(y - y.mean()))

sxx = np.sum((x - x.mean())**2)

syy = np.sum((y - y.mean())**2)

a = sxy / sxx

b = y.mean() - a*x.mean()

r = sxy / np.sqrt(sxx * syy)

n = len(num_stores)

ESS = n * sxy**2 / sxx

TSS = n * syy

R2 = ESS / TSS

max_population = 15000000

max_stores = 8000

plt.xlim(0, max_population)

plt.ylim(0, max_stores)

plt.xlabel('population')

plt.ylabel('number of stores')

plt.scatter(population, num_stores, label='#stores')

x = np.linspace(0, max_population)

plt.plot(x, a*x + b, color='m')

plt.text(1000000, 7600, 'a = ' + str(a))

plt.text(1000000, 7300, 'b = ' + str(b))

plt.text(1000000, 7000, 'r = '+str(r))

plt.text(1000000, 6700, 'R2 = '+str(R2))

plt.show()

人口との関係

コンビニ店舗数と人口の散布図と回帰式を以下に示す。

相関係数が極めて高いのは当然で、やはりコンビニの出店計画には人口ファクターが強く作用していることがわかる。

定数項bが店舗数のオーダーに比べてほぼゼロというのも興味深い。

係数aの値からは5万人弱で1店舗ということになるが、人口10万くらいの都市で2店舗しかないことになり、ちょっと少ないような気がする。もしかすると、一定規模以上の市町村単位や都市単位くらいで層化して出店計画を立てているのかもしれない。

面積との関係

店舗数と面積の関係

次に店舗数と都道府県面積の関係を見てみた。

この結果はかなりはずれで、データを見ても人口が少なく面積が群を抜いて大きい北海道の影響を大きく受けている。

ここで、面積が極端に大きい北海道（83423.83㎢～1位、2906店舗～5位）と面積が小さいが集積度が極端に高い東京都（2193.96㎢～45位、7003店舗～1位）の2つを除いて計算してみる。

これはさらにおかしな結果で、面積が小さいほど店舗数が多いことになる。

考えてみれば、集客数を期待するなら人口が集積している地域が有利だから、人口密度に比例する可能性を考えた方がいいのかもしれない。もし面積が小さい県の方が集積度が高いと想定すると、面積だけを取り出したときに逆の関係になるとも考えられるが、相関係数や決定係数が小さすぎるので考察は難しい。

人口密度

以下は人口密度との関係。

今度はかなりきれいに相関の高さが出ている。

直接的な計算式に入れているかどうかわからないが、GISなどで出店計画を立てるとしたら、人口密度の高いエリアを選んでいくだろうことが想定される。

ただ、店舗数は人口などといった売り上げに直結するデータから導かれるのが普通で、人口密度が高くても人口が少なければ出店インセンティブにはならない。

人口と人口密度の関係

試しに人口を説明変数、人口密度を被説明変数として両者の関係を見てみると、驚くことに「人口が多いほど人口密度が高くなる（あるいはその逆）」という関係になる。

ここから先は人口論や地域論になりそうなので置いておくが、少なくとも日本においては、「狭いところほど人が集まっている傾向がある」ということになりそうである。

もちろんこれは他の国でも一般に当てはまることかもしれないが、朝のラッシュ時に特定の車両に無理やり乗り込んでいる割に離れた車両がすいているとか、1本電車を遅らせたらガラガラだったとか、そのあたりの行動パターンを見ていると、何となく日本に特有のような気がする。

道路延長との関係

最後に道路延長との関係を見てみる。

東京のように稠密な都市は例外とすると、概ね関係はありそうである。ただし相関係数、決定係数は高くない。

コンビニ店舗が道路の利便性に依っていることは推測できるが、やはり人口という売り上げ直結のデータに比べると関係は弱い。

高速道路の延長についても見てみたが、こちらはほとんど関係は見られなかった。

ただ、高速道路の延伸に伴ってコンビニエンスストアの店舗数が伸びているようであり、マクロな延長というよりも物流上のインパクトが大きいことは予想される。

相関係数

2019-09-29 / tau / コメントする

相関係数の定義

相関係数は、多数のデータの組がどの程度線形に近い性質を持つかを表す値で、以下で定義される。

(1) $\begin{eqnarray*} r &=& \frac{{\rm Cov}(X, Y)}{\sqrt{V(X) \cdot V(Y)}} \\ &=& \frac{E(X - \overline{X})(Y - \overline{Y})} {\sqrt{ ( E\left[ (X - \overline{X})^2 \right] E\left[ (Y - \overline{Y})^2 \right] }} \\ &=& \frac{\displaystyle \sum_{i=1}^{n}(x_i - \overline{x})(y_i - \overline{y})} {\left[ \displaystyle \sum_{i=1}^{n} (x_i - \overline{x})^2 \displaystyle \sum_{i=1}^{n} (y_i - \overline{y})^2 \right]^{1/2}} \end{eqnarray*}$

相関係数の線形変換

変数を線形変換した場合の相関係数は分散・共分散の性質から、以下のようになり、元のままか符号が反転する。

(2) $\begin{eqnarray*} r' &=& \frac{{\rm Cov}(aX+b, cY+d)}{\sqrt{V(aX+b) \cdot V(cY+d)}} \\ &=& \frac{ac{\rm Cov}(X, Y)}{\sqrt{a^2 V(X) \cdot c^2 V(Y)}} \\ &=& \frac{ac}{|ac|} \cdot r \end{eqnarray*}$

完全線形関係の場合の相関係数

XとYが完全な線形関係にある場合、相関係数は1または-1になる。このとき、傾きの大きさや、平行移動は相関係数に影響しない。

(3) $\begin{eqnarray*} r &=& \frac{ {\rm Cov}(X, aX+b) } { \left[ V(X) \cdot V(aX+b) \right]^{1/2} } \\ &=& \frac{ a {\rm Cov}(X, X) } { \left[ V(X) \cdot a^2 V(X) \right]^{1/2} } \\ &=& \frac{a}{\sqrt{a^2}} \cdot \frac{ {\rm Cov}(X, X) } { \left[ V(X) \cdot V(X) \right]^{1/2} } \\ &=& \frac{a}{|a|} \cdot \frac{ V(X) } { V(X) } \\ &=& 1 \; {\rm or} \; -1 \end{eqnarray*}$

いろいろな分布の相関係数

完全な線形関係

以下のコードで確認。

相関係数は傾きや平行移動に対して影響を受けず、増加関数なら1、減少関数なら－1になることがわかる。

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

def cor(x, y):
    return np.sum((x - x.mean())*(y - y.mean())) / len(x) / x.std() / y.std()

x = np.linspace(0.2, 0.8, num=21)

y01 = np.array([x for x in x])
y02 = np.array([x * 0.5 for x in x])
y03 = np.array([x + 0.2 for x in x])
y04 = np.array([1 - x for x in x])

plt.axes().set_aspect('equal')
plt.xlim(0, 1)
plt.ylim(0, 1)

plt.scatter(x, y01, label='cor=' + str(cor(x, y01)))
plt.scatter(x, y02, label='cor=' + str(cor(x, y02)))
plt.scatter(x, y03, label='cor=' + str(cor(x, y03)))
plt.scatter(x, y04, label='cor=' + str(cor(x, y04)))

plt.legend()

plt.show()

import numpy as np

import numpy.random as rnd

import matplotlib.pyplot as plt

def cor(x, y):

return np.sum((x - x.mean())*(y - y.mean())) / len(x) / x.std() / y.std()

x = np.linspace(0.2, 0.8, num=21)

y01 = np.array([x for x in x])

y02 = np.array([x * 0.5 for x in x])

y03 = np.array([x + 0.2 for x in x])

y04 = np.array([1 - x for x in x])

plt.axes().set_aspect('equal')

plt.xlim(0, 1)

plt.ylim(0, 1)

plt.scatter(x, y01, label='cor=' + str(cor(x, y01)))

plt.scatter(x, y02, label='cor=' + str(cor(x, y02)))

plt.scatter(x, y03, label='cor=' + str(cor(x, y03)))

plt.scatter(x, y04, label='cor=' + str(cor(x, y04)))

plt.legend()

plt.show()

線形性が強い関係

線形関数の値に対して、乱数でばらつきを与えた場合の相関係数の違いを示す。ばらつきが大きい方が相関係数は小さくなる。

y11 = np.array([x for x in x])
y12 = np.array([x + rnd.rand() * 0.2 * rnd.choice([-1, 1]) for x in x])
y13 = np.array([x + rnd.rand() * 0.4 * rnd.choice([-1, 1]) for x in x])

y11 = np.array([x for x in x])

y12 = np.array([x + rnd.rand() * 0.2 * rnd.choice([-1, 1]) for x in x])

y13 = np.array([x + rnd.rand() * 0.4 * rnd.choice([-1, 1]) for x in x])

負の線形性が強い場合は、相関係数がマイナスになる。

y21 = np.array([1 - x for x in x])
y22 = np.array([1 - x + rnd.rand()*0.2*rnd.choice([-1, 1]) for x in x])
y23 = np.array([1 - x + rnd.rand()*0.4*rnd.choice([-1, 1]) for x in x])

y21 = np.array([1 - x for x in x])

y22 = np.array([1 - x + rnd.rand()*0.2*rnd.choice([-1, 1]) for x in x])

y23 = np.array([1 - x + rnd.rand()*0.4*rnd.choice([-1, 1]) for x in x])

放物線（強い関係があるのに相関が低いケース）

以下のような放物線では、XとYにきちんとした数学的関係があるのに、相関係数がゼロに近くなる。

相関係数はXとYが単調増加／単調減少の度合いが強いほど、また線形関係に近いほど1に近くなるが、それ以外の関係が強い場合にはそれを補足できない場合がある。

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

def cor(x, y):
    return np.sum((x - x.mean())*(y - y.mean())) / len(x) / x.std() / y.std()

x = np.linspace(0, 1.0, num=21)

y1 = np.array([4 * (x - 0.5)**2 for x in x])

plt.axes().set_aspect('equal')
plt.xlim(0, 1)
plt.ylim(0, 1)

plt.scatter(x, y1)

plt.text(0.05, 0.95, 'cor=' + str(cor(x, y1)))

plt.show()

import numpy as np

import numpy.random as rnd

import matplotlib.pyplot as plt

def cor(x, y):

return np.sum((x - x.mean())*(y - y.mean())) / len(x) / x.std() / y.std()

x = np.linspace(0, 1.0, num=21)

y1 = np.array([4 * (x - 0.5)**2 for x in x])

plt.axes().set_aspect('equal')

plt.xlim(0, 1)

plt.ylim(0, 1)

plt.scatter(x, y1)

plt.text(0.05, 0.95, 'cor=' + str(cor(x, y1)))

plt.show()

反比例（負の線形性に見えてしまう場合）

以下は反比例関数の場合。

関数の形状や範囲に寄るが、この場合は相関係数の絶対値が0.8以上と1に近く、これだけ見ると負の線形性が強そうに見える。

x = np.linspace(0.1, 1.0, num=21)
y2 = np.array([0.1 / x for x in x])

1 2	x = np.linspace(0.1, 1.0, num=21) y2 = np.array([0.1 / x for x in x])

対数関数（正の線形性に見えてしまう場合）

対数関数の場合。この場合は0.9以上とかなり強い線形性を示唆している。

x = np.linspace(0.1, 1.0, num=21)
y3 = np.array([np.log(x*20-1)/4 for x in x])

1 2	x = np.linspace(0.1, 1.0, num=21) y3 = np.array([np.log(x*20-1)/4 for x in x])

相関係数に関する注意

本来の関係との乖離

先にみたように、線形関係ではないが数学的な関係を持つ場合に、相関係数からは全く関係がない、元の関係とは異なり線形関係を持つ、といった解釈になることがある。

相関係数が高い場合に、線形回帰式などで物事を予測する際には注意が必要。

変数間に解析的な関係が見いだせるならそれを重視すべきであり、よしんばそれがわからないにしても、定義域の範囲で「ある程度は当たる」程度に考えておくべきか。

因果関係

堂々と間違えられるケースが、「科学的な」記事やマスメディアなどでよくみられる。

相関係数は「変数間の単調な増加／減少傾向が強いかどうか」だけを示すもので、必ずしも因果関係を示唆しない。

気温が高いとビールはよく売れるがおでんは売れない。その二つに負の相関があるからといって、「ビール好きはおでんが嫌い」と言えないが、形を変えてこのような解釈がなされる恐れがある。

もともとのメカニズムで因果関係が示唆されていて、その上で相関係数の大きさを論じるなら意義もあるが、その場合でも、事象に対する寄与度などをよく考えておかないと「それだけが原因」と考えるような間違いを犯すことになる。

回帰分析～単回帰

2019-09-09 / tau / コメントする

概要

$(X,\: Y) = (x_i,\: y_i) \; (i=1 \ldots n)$ のデータからなる複数の標本に対して、線形式で最も説明性の高いものを求める。

そのために、各データの $x_i$ に線形式を適用して得られた値 $\hat{y_i} = a x_i + b$ と $y_i$ の差（残差）の平方和が最小となるような係数 $a,\: b$ を求める（最小二乗法）。

定式化

(1) $\begin{equation*} \sum_{i=1}^n (\hat{y_i} - y_i)^2 \; \rightarrow \; \min \quad {\rm where} \; \hat{y_i} = a x_i + b \end{equation*}$

ここで残差を最小化する $a,\: b$ を求めるために、それぞれで偏微分する。

(2) $\begin{equation*} \left\{ \begin{array}{l} \displaystyle \frac{\partial}{\partial a} \sum_{i=1}^n (a x_i + b - y_i)^2 = 0 \\ \displaystyle \frac{\partial}{\partial b} \sum_{i=1}^n (a x_i + b - y_i)^2 = 0 \end{array} \right. \quad \Rightarrow \quad \left\{ \begin{array}{l} \displaystyle \sum_{i=1}^n 2(a x_i + b - y_i)x_i = 0 \\ \displaystyle \sum_{i=1}^n 2(a x_i + b - y_i) = 0 \end{array} \right. \end{equation*}$

展開して

(3) $\begin{equation*} \left\{ \begin{array}{l} \displaystyle a \sum_{i=1}^n x_i^2 + b \sum_{i=1}^n x_i - \sum_{i=1}^n x_i y_i = 0 \\ \displaystyle a \sum_{i=1}^n x_i + bn - \sum_{i=1}^n y_i = 0 \end{array} \right. \end{equation*}$

各和の記号を新たに定義し、行列表示で表すと、

(4) $\begin{equation*} \left[ \begin{array}{cc} S_{xx} & S_x \\ S_x & n \end{array} \right] \left[ \begin{array}{c} a \\ b \end{array} \right] = \left[ \begin{array}{c} S_{xy} \\ S_y \end{array} \right] \end{equation*}$

回帰式の係数

式(4)を解くと以下を得る。

(5) $\begin{equation*} \left[ \begin{array}{c} a \\ b \end{array} \right] = \left[ \begin{array}{c} \dfrac{n S_{xy} - S_x S_y}{n S_{xx} - {S_x}^2} \\ \dfrac{S_{xx} S_y - S_x S_{xy}}{n S_{xx} - {S_x}^2} \end{array} \right] \end{equation*}$

あるいは式(5)については以下のようにも表せる。ただし以下の式で、 $\sigma_{XY} = {\rm Cov}(X, Y)$ 、 ${\sigma_X}^2 = {\rm Var}(X)$

(6) $\begin{equation*} a = \frac{\displaystyle \sum_{i=1}^n (x_i - \overline{X})(y_i - \overline{Y})}{\displaystyle \sum_{i=1}^n (x_i - \overline{X})^2} =\frac{\sigma_{XY}}{{\sigma _X}^2} \end{equation*}$

(7) $\begin{equation*} b = \overline{Y} - a \overline{X} \end{equation*}$

決定係数

全変動・回帰変動・残差変動

$x_i$ の回帰式による推定値を $\hat{y}_i = a x_i + b$ で表す。

ここで、 $y_i,\: \hat{y},\: \overline{Y}$ によって、以下の変動が定義される。

$y_i - \overline{Y}$: 全変動。平均からのデータのばらつき。
$\hat{y}_i - \overline{Y}$: 回帰変動。平均からのばらつきのうち、回帰式によって説明される部分。
$y_i - \hat{y}_i$: 残差変動。平均からのばらつきのうち、回帰式によって説明されない部分。

当然、[全変動]＝[回帰変動]＋[残差変動]となることはすぐに確認できる。

また、これら3つの変動の2乗和をTSS(Total Sum of Squares)、ESS(Explained Sum of Squares、RSS(Residual Sum of Squares)として、以下のように定義する。

(8) $\begin{eqnarray*} TSS &=& \sum_{i=1}^n (y_i - \overline{Y})^2 \\ ESS &=& \sum_{i=1}^n (\hat{y}_i - \overline{Y})^2 \\ RSS &=& \sum_{i=1}^n (y_i - \hat{y}_i)^2 \end{eqnarray*}$

このとき、それぞれについて以下のように表せる。

(9) $\begin{eqnarray*} TSS &=& n {\sigma _Y}^2 \\ ESS &=& n \frac{{\sigma_{XY}}^2}{{\sigma_x}^2} \\ RSS &=& n {\sigma_Y}^2 - n \frac{{\sigma_{XY}}^2}{{\sigma_X}^2} \end{eqnarray*}$

これらから、全変動・回帰変動・残差変動の2乗和についても、以下の関係が成り立つことがわかる。

(10) $\begin{equation*} TSS = ESS + RSS \end{equation*}$

補足～ESS

式(9) のESSの確認。以下のように、線形変換された分散から導ける。

(11) $\begin{eqnarray*} ESS &=& \sum_{i=1}^n \left( a x_i +b - (a \overline{X} + b) \right)^2 \\ &=& n {\rm Var}(a X + b) \\ &=& na^2{\rm Var}(X) \\ &=& n a^2 {\sigma_x}^2 \\ &=& n \left( \frac{\sigma_{XY}}{{\sigma_x}^2} \right)^2 {\sigma_x}^2 \\ &=& n \frac{{\sigma_{XY}}^2}{{\sigma_x}^2} \end{eqnarray*}$

あるいは、以下のように地道に計算しても確認できる。

(12) $\begin{eqnarray*} \sum_{i=1}^n (\hat{y}_i - \overline{Y}) &=& \sum_{i=1}^n ({\hat{y}_i}^2 - 2 \hat{y}_i \overlne{Y} + \overline{Y}^2) \\ &=& \sum_{i=1}^n (\hat{y}_i^2 - 2(a x_i + b) \overline{Y} +\overline{Y}^2) \\ &=& \sum_{i=1}^n (\hat{y}_i^2 - 2a x_i \overline{Y} - 2b\overline{Y} + \overline{Y}^2) \\ &=& \sum_{i=1}^n \hat{y}_i^2 - 2an \overline{X} \; \overline{Y} -2nb \overline{Y} + n \overline{Y}\end{eqnarray*} \\$

ここで

(13) $\begin{eqnarray*} \sum_{i=1}^n \hat{y}_i^2 &=& \sum_{i=1}^n (ax_i + b)^2 \\ &=&\sum_{i=1}^n (a^2 x_i^2 + 2ab x_i + b^2) \\ &=& a^2 \sum_{i=1}^n x_i^2 + 2abn \overline{X} + nb^2 \end{eqnarray*}$

これを代入して、

(14) $\begin{eqnarray*} ESS &=& a^2 \sum_{i=1}^n x_i^2 + 2nab \overline{X} + nb^2 - 2na \overline{X} \; \overline{Y} - 2nb \overline{Y} + n \onerline{Y} \end{eqnarray*}$

さらに $\sum_{i=1}^n x_i^2 = \n \sigma_X^2 + n \overline{X}^2$ を考慮して、

(15) $\begin{eqnarray*} ESS &=& n a^2 {\sigma_X}^2 + n a^2 \overline{X}^2 + 2nab \overline{X} + nb^2 -2na\overline{X}\;\overline{Y} - 2nb\overline{Y} + n\overline{Y}^2 \\ &=& n a^2 {\sigma_X}^2 + n(a \overline{X} - \overline{Y})^2 + 2nab\overline{X} + nb^2 - 2nb\overline{Y} \\ &=& n a^2 {\sigma_X}^2 + nb^2 + 2nab\overline{X} + nb^2 - 2nb\overline{Y} \\ &=& n a^2 {\sigma_X}^2 + 2nb(b + a\overline{X} - \overline{Y}) \\ &=& n a^2 {\sigma_X}^2 \\ &=& n \left( \frac{\sigma_{XY}}{{\sigma_X}^2} \right)^2 {\sigma_X}^2 \\ &=& n \frac{{\sigma_{XY}}^2}{{\sigma_X}^2} \end{eqnarray*}$

補足～RSS

式(9) のRSSの確認。

(16) $\begin{eqnarray*} RSS &=& \sum_{i=1}^n (y_i - \hat{y}_i)^2 = \sum_{i=1}^n (y_i - a x_i - b)^2 \\ &=& \sum_{i=1}^n (y_i - a x_i - \overline{Y} + a \overline{X})^2 \\ &=& \sum_{i=1}^n \left(y _i - \overline{Y} - a (x_i - \overline{X}) \right) ^2 \\ &=& \sum_{i=1}^n \left(y _i - \overline{Y} - \frac{\sigma_{XY}}{{\sigma_X}^2} (x_i - \overline{X}) \right) ^2 \\ &=& n {\sigma_Y}^2 - 2n \frac{\sigma_{XY}}{{\sigma_X}^2} \sigma_{XY} + n \frac{{\sigma_{XY}}^2}{{\sigma_X}^4} {\sigma_X}^2 \\ &=& n {\sigma_Y}^2 - 2n \frac{{\sigma_{XY}}^2}{{\sigma_X}^2} + n \frac{{\sigma_{XY}}^2}{{\sigma_X}^2} \\ &=& n {\sigma_Y}^2 - n \frac{{\sigma_{XY}}^2}{{\sigma_X}^2} \end{eqnarray*}$

決定係数

決定係数は、全変動のうち回帰変動でどれだけ説明できるか（残差変動が少ないか）を表したもので、通常 $R^2$ で表される。

(17) $\begin{align*} R^2 = \frac{ESS}{TSS} = \frac{TSS - RSS}{TSS} = 1 - \frac{RSS}{TSS} \end{align*}$

この値は、残差変動が小さいほど1に近づき、 $RSS = TSS$ すなわち全変動が回帰変動で全く説明できないときに0となる。

ロシアンルーレットの確率

2019-07-04 / tau / コメントする

基本

ロシアンルーレットは、賭の順番問題や条件付き確率の問題としてよく見かける。まず、一番基本的なルールで考える。

1人の相手とロシアンルーレットで命をかけることになった。6発入りの拳銃に1発だけ弾を込め、最初にシリンダーを高速で回転させる。シリンダーが止まったところで、相手と交互に引き金を引いていくとして、まず自分から始めるのと相手に譲るのと、どちらが生き残る可能性が高いか。

n回目に実弾に当たる確率を $P_n \; (n = 1~6)$ とすると、

1発目で当たる確率は、6つの弾倉のどこかに実弾が入っているので、

(1) $\begin{equation*} P_1 = \frac{1}{6} \end{equation*}$

2発目で当たる確率は、1発目で実弾に当たらず、2発目は残り5つの弾倉のどこかに実弾が入っているので、

(2) $\begin{equation*} P_2 = \frac{5}{6} \times \frac{1}{5} = \frac{1}{6} \end{equation*}$

以下、6発目まで同じように考えて、

(3) $\begin{eqnarray*} P_3 &=& \frac{5}{6} \times \frac{4}{5} \times \frac{1}{4} = \frac{1}{6} \\ P_4 &=& \frac{5}{6} \times \frac{4}{5} \times \frac{3}{4} \times \frac{1}{3} = \frac{1}{6} \\ P_5 &=& \frac{5}{6} \times \frac{4}{5} \times \frac{3}{4} \times \frac{2}{3} \times \frac{1}{2} = \frac{1}{6} \\ P_6 &=& \frac{5}{6} \times \frac{4}{5} \times \frac{3}{4} \times \frac{2}{3} \times \frac{1}{2} \times \frac{1}{1}= \frac{1}{6} \end{eqnarray*}$

すなわち、通常のロシアンルーレットでは、1～6回目のどこで実弾に当たる確率も等しく1/6で、順番の後先に有利性はない。

毎回シリンダーを回転させる場合

1発だけ弾を込めた拳銃で、2人でロシアンルーレットで命をかける。ただし、それぞれが自分の番になる毎にシリンダーを回転させて、ランダムに止まったところで引き金を引く。自分が最初に打つのと、相手から打たせるのと、どちらを選ぶか。

毎回シリンダーを回転させるので、そこまで発砲されていなければ、実弾である確率は毎回1/6。 $n$ 回目に実弾に当たる確率を $p_n$ とすると、 $n-1$ 回は当たらないことを考慮して、

(4) $\begin{equation*} p_n = \left( \frac{5}{6} \right) ^{n-1} \times \frac{1}{6} \end{equation*}$

まず、この確率はnが大きくなるほどゼロに近づく。たとえば1回めで当たる確率は約16.7%、10回目で当たる確率は約3.2%、20回目で当たる確率は0.5%、40回目だと約0.01%程度となる。

次に、 $n$ 回目までに当たる確率を $P_n$ とすると、これは各回の確率を累積していくことで求められるので、等比数列の和より、

(5) $\begin{equation*} P_n = \sum_{i=1}^{n} \left( \frac{5}{6} \right) ^{n-1} \times \frac{1}{6} = \frac{1}{6} \frac{1-\left( \dfrac{5}{6} \right)^n}{1 - \dfrac{5}{6}} = 1-\left( \dfrac{5}{6} \right)^n \end{equation*}$

この累積確率は $n \rightarrow \infty$ のとき1に近づくが、10回目で約83.8%、20回目で約97.4%、30回目で約99.6%、40回目で約99.9%となり、100回も打ってその時点で生き延びている確率はほとんどゼロに近い。

次に交互に打っていく場合を考える。

先攻で打つ方になった場合、1回目、3回目・・・と奇数回に打つことになるので、その累積確率を $P_{1 \cdot n}$ とすると、

(6) $\begin{eqnarray*} P_{1\cdot 2n-1} &=& \sum_{k=1}^{n} \left( \frac{5}{6} \right)^{2(k-1)} \times \frac{1}{6} = \sum_{k=1}^{n} \frac{1}{6} \left( \frac{25}{36} \right)^{k-1} = \frac{1}{6} \cdot \frac{1 - \left( \dfrac{25}{36}\right) ^n }{1 - \dfrac{25}{36}} \\ &=& \frac{6}{11} \left( 1 - \left( \dfrac{25}{36}\right) ^n \right) \end{eqnarray*}$

後攻で打つ方になった場合は、2回目、4回目・・・と偶数回に打つことになるので、その累積確率を $P_{2 \cdot n}$ とすると、

(7) $\begin{eqnarray*} P_{2\cdot 2n} &=& \sum_{k=1}^{n} \left( \frac{5}{6} \right)^{2k-1} \times \frac{1}{6} = \sum_{k=1}^{n} \frac{5}{36} \left( \frac{25}{36} \right)^{k-1} = \frac{5}{36} \cdot \frac{1 - \left( \dfrac{25}{36}\right) ^n }{1 - \dfrac{25}{36}} \\ &=& \frac{5}{11} \left( 1 - \left( \dfrac{25}{36}\right) ^n \right) \end{eqnarray*}$

同じ $n$ に対して常に $P_{1 \cdot n} / P_{2 \cdot n} = 1.2$ となり、先に打つのを選んだほうが実弾に当たる確率は高く、この比率は $n \rightarrow \infty$ としたときの結果と同じ。

したがって、「お先にどうぞ」と相手に最初を譲ったほうが、少しでも生き残る確率が高くなる。また、3人以上で順番に打っていく場合は、できるだけ後の順番で打つ方が実弾に当たる確率が低くなる。

条件付きの問題

2発の弾を並べて込める

ロシアンルーレットと少し違うが、条件付きの問題。

筋悪の金貸しから借りたかなりの借金を返せなくなり、命を取られようとしているが、「チャンスをやる」と言われた。彼は6発入りの拳銃の弾倉に2 発の弾を並べて込め、シリンダーを回す。シリンダーがランダムな位置で止まったところで、彼は言う。「今から空に向かって私が1回引き金を引く。その後にお前のこめかみに銃口を当てて打ってもらうが、そのまま打ってもいいし、もう一度シリンダーを回して止まった位置で打ってもいい」

1発目が空砲の場合

金貸しが銃を空に向けて引き金を引くと空砲だった。そのまま打つか、シリンダーを回転させるか、どちらを選択すべきか。

まず、シリンダーをもう一度回転させてから打つ場合は、6弾倉中2発残っているので、当たる確率は1/3。逆に助かる率は2/3となる。

一方シリンダーを回転させない場合は、手渡された拳銃の弾の状態を、円形に並んだ弾倉を一列で表し、下表のように整理する。ただし、左欄の●は実弾がある場所、○は実弾が入っていない場所で、1番目は金貸しが打ったときの位置、2番目は手渡された時の位置。右欄は2発目が実弾の場合は●、空砲の場合は○とする。

●●○○○○	対象外
○●●○○○	●
○○●●○○	○
○○○●●○	○
○○○○●●	○
●○○○○●	対象外

金貸しが打ったときは空砲だったので、1番目が実弾の事象は対象外。このとき、実弾に当たらない確率は3/4となる。

すなわち、もう一度シリンダーを回転させるよりも、そのまま続けて打った方が助かる可能性が高い。

1発目が実弾の場合

もし、金貸しが最初に空へ向けて打ったときに実弾だった場合はどうか？

この場合に続けて打つと、1発目が実弾である2ケースに対して、空砲で助かる率は1/2。

シリンダーを回転させてから打つと、実弾は残り1発なので、当たらない確率は5/6で、回転させて打った方が生存確率は高くなる。

条件に応じた選択と確率

以上の結果を整理すると下表の通りとなり、1発目が空砲の場合と実弾の場合で、その後の選択による生存確率の高い方が逆転している。

空砲も実弾もその位置が連続しているため、1発目が空砲の場合は続けて空砲の確率が高く、1発目が実弾の場合は次も実弾の確率が高そうだが、1発目が実弾の場合は弾が1つ消費されるので、その影響で結果が逆転している。

1発目	選択	生存確率
空砲	シリンダー回転	1/3 (0.33)
空砲	続けて打つ	3/4 (0.75)
実弾	シリンダー回転	5/6 (0.83)
実弾	続けて打つ	1/2 (0.50)

実弾2発をランダムに込める

上記の金貸しの問題では実弾2発を並べて込めたが、これをランダムな位置に2発込めた場合はどうなるか。

この場合、2発の実弾の装填状況は下表のようになる。

1発目が→	空砲	実弾
●●○○○○	対象外	●
●○●○○○	対象外	○
○●●○○○	●	対象外
●○○●○○	対象外	○
○●○●○○	●	対象外
○○●●○○	○	対象外
●○○○●○	対象外	○
○●○○●○	●	対象外
○○●○●○	○	対象外
○○○●●○	○	対象外
●○○○○●	対象外	○
○●○○○●	●	対象外
○○●○○●	○	対象外
○○○●○●	○	対象外
○○○○●●	○	対象外

1発目が空砲だった場合

シリンダーを回転させた場合の生存確率は、4/6 = 2/3 (0.667)。

続けて打つ場合は、1発目が実弾のケースを除いて2発目が空砲の場合なので、6/10 = 3/5 (0.6)。

1発目が実弾だった場合

シリンダーを回転させた場合の生存確率は、5/6 (0.833)。

続けて打つ場合は、1発目が空砲のケースを除いて2発目が空砲の場合なので、4/5 (0.8)。

1発目が空砲でも実弾でも、シリンダーを回転させた方が生き残る確率が高くなる。

条件付確率～男女問題

2019-07-01 / tau / 4件のコメント

概要

条件付き確率の代表的な問題として、2人の子の性別問題をよく見かけるが、前提条件によって答えが異なることを明確にしている例は少ない。

また、この問題がベン図を用いた分析や、ベイズの定理の親和性についても考えてみる。

兄妹の場合

まず準備として、以下の問題を考える。

2人兄弟の上の子が男であった場合、下の子が女である確率はいくらか。

普通に考えれば、男か女どちらかなのだから、確率は1/2と考えるが、まずは生真面目にすべての事象を並べてみる。

conditional-probability-brother-cases

上の子が兄か姉かに関係なく下の子が女である確率は1/2と直感的に考えられるが、それが確認できた。

一応、以下の様なベン図を考えてみる。図中、色付けしたところが条件下での事象で、濃い色が事後確率として求めたい部分。

conditional-probability-brother-venn-diagram

ここで記号E/Yはelder/younger、M/Fはmail/femailを表す。この場合、上の子が男（事象EM）のもとでの下の子が女（事象EMかつYF）ということになり、EMの要素が(男, 男)、(男, 女)の2つに対してEMかつYFの要素は(男,女)だけなので、確率は1/2となる。

なお、ベン図の代わりに場合分けの表で整理してみると、以下の様になる。
conditional-probability-brother-table

これを敢えてベイズの定理で書くと以下の様になるが、まあ当たり前という感じ。

(1) $\begin{equation*} P(YF|EM) = \frac{P(EM \cap YF)}{P(EM)} = \frac{1}{2} \end{equation*}$

これをベイズの定理の別の表現で書くと以下の様になる。

(2) $\begin{equation*} P(YF|EM) = \frac{P(EM|YF)P(YF)}{P(EM)} = \frac{\dfrac{1}{2} \cdot \dfrac{1}{2}}{\dfrac{1}{2}}=\frac{1}{2} \end{equation*}$

ただし $P(EM|YF)=1/2$ としているが、この項の考え方は結局求めたい確率と同じ道筋で考えなければならない。

いずれにしてもこの問題は、第一子が男であろうが女であろうが、第二子が男か女かと問われれば、確率1/2でどちらか、という自然な結果となり、ベン図や場合分け表、ベイズの定理を使ったとしても、あまり恩恵は得られない。

区別なしの場合

次に以下の様な問題を考える。

2人の子のうち、1人が男の場合、もう1人が女である確率はいくらか。

これも普通に考えれば、男か女どちらかなのだから、確率は1/2と考えてしまいそうになるが、まじめに事象を考えてみる。

conditional-probability-2children-cases

この場合は1人が男とわかっているので、(姉, 妹)の事象は存在しなくなり、もう1人が女というだけで姉か妹かの区別はないことから、確率は2/3となる。

これをベン図で描いてみると以下のようになるが、男と女の組み合わせが2通りあることを確認する点では、上のような整理のほうがわかりやすい。

conditional-probability-2children-venn-diagram-1

場合分けの表は以下の通りになる。

conditional-probability-2children-table-1

これをベイズの定理で表現すると以下のようになる。

(3) $\begin{eqnarray*} P(F|M) &=& \frac{P(F \cap M)}{P(M)} = \frac{P(FM) + P(MF)}{P(MM) + P(MF) + P(FM)} \\ &=& \frac{\dfrac{1}{4} + \dfrac{1}{4}}{\dfrac{1}{4} + \dfrac{1}{4} + \dfrac{1}{4}} = \frac{2}{3} \end{eqnarray*}$

条件付き確率から求める方法は、式(2)と同じ理由でメリットがない。

順番がある場合

次に以下の様な問題を考える。

2人の子がいる家庭を訪ねた。最初に挨拶に出てきたのが男の子だった場合、もう1人が女の子である確率はいくらか。

この場合は上の子と下の子を区別していないので、一つ前の問題と同じ答えになりそうだが、事象を列挙してみる。

conditional-probability-2children-by-order-cases

この問題でも、2人とも女のケースが除かれるが、兄弟構成に対して挨拶の順番がそれぞれ2通りずつあり、2人とも男の場合は兄、弟のいずれも最初に挨拶しうるので前提条件の事象に2つ含まれる。結果として、2人目が女の子である確率は1/2となる。

このように、2人兄弟の性別に関する問題でも、条件の付し方によって答えが異なってくるが、いずれの問題もベイズの定理を駆使するというものではなさそうである。

3人兄弟の場合

問題を以下の様に拡張する。

3人の子のうち2人が男の子の場合、残る1人が女の子である確率はいくらか。

これについても以下のように事象を列挙し、確率3/4を得る。

conditional-probability-3children-cases

ベイズの定理で表すと以下のようになるが、この場合もまじめに事象のパターンを列挙した考えるほうが素直。

(4) $\begin{eqnarray*} P(F|MM) &=& \frac{P(F \cap MM)}{P(MM)} \\ &=& \frac{P(MMF) + P(MFM) + P(FMM)}{P(MMM) + P(MMF) + P(MFM) + P(FMM)} \\ &=& \frac{\dfrac{1}{8} + \dfrac{1}{8} + \dfrac{1}{8}}{\dfrac{1}{8} + \dfrac{1}{8} + \dfrac{1}{8} + \dfrac{1}{8}} = \frac{3}{4} \end{eqnarray*}$

ベイズの定理～難病検査

2019-07-01 / tau / コメントする

概要

ベイズの定理で典型例として説明される問題。ある病気の検査で陽性が出た場合に罹患している確率を求める。類似の問題として、鉱脈を探るのに試掘結果から鉱脈が存在する確率を求める問題がある。

ロジックとしては、「誤判定の可能性がある調査・検査の結果から、目的とする事実の確率を求める」という構図で共通している。

難病問題

病気の検査の問題を例示する。

ある難病の検査結果が陽性だったとき、その病気に罹患している確率を求めよ。ただし以下の情報が与えられている。

その難病は1000人に1人が罹患している
その難病にかかっている場合、90%の確率で検査結果は陽性になる
その難病にかかっていなくても、20%の確率で検査結果は陽性になる

ここで以下の記号を定義する。

全体としての難病の罹患率→ $P(D)=0.001$
難病にかかっている場合に陽性になる確率→ $P(P|D)=0.9$
難病にかかっていない場合に陽性になる確率→ $P(P| \overline{D})=0.2$

このとき、ベイズの定理により以下の確率を得る。

(1) $\begin{eqnarray*} P(D|P) &=& \frac{P(P|D) P(D)}{P(P|D) P(D) + P(P|\overline{D}) P(\overline{D})} \\ &=& \frac{0.9 \times 0.001}{0.8 \times 0.001 + 0.2 \times 0.999} \\ &\approx& 0.004484 \end{eqnarray*}$

検査結果が陽性であっても、この難病にかかっている確率は0.5%未満ということになる。検査しなければ1000人に1人の罹患率なので0.1%なので、検査によってその確からしさが約4.5倍となるが、事前確率が小さい場合、検査の精度を上げても劇的な判定効果は望めない。

なお、この約4.5倍という値は $P(P|D) / P(P|\overline{D}) = 4.5$ に相当し、事前確率 $P(D)$ の値が小さいほど4.5に近づく(ベイズの定理の解釈参照)。

なお、2020年、世界的に大きな影響を及ぼしたCOVID-19（新型コロナウイルス）のPCR検査を題材にした記事をこちらにまとめている。

ベイズの定理の解釈

2019-07-01 / tau / コメントする

ベイズの定理の変形

ベイズの定理は、一般に以下のように表現される。

(1) $\begin{equation*} P(T_k |C) = \frac{P(C|T_k)P(T_k)}{\displaystyle \sum_{i=1}^n P(C|T_i)P(T_i)} \end{equation*}$

または $T$ と $\overline{T}$ の2事象の場合は、

(2) $\begin{equation*} P(T|C) = \frac{P(C|T)P(T)}{P(C|T)P(T) + P(C|\overline{T})P(\overline{T})} \end{equation*}$

ここで次のように確率の記号を定義する。また、イメージを掴みやすくするために、難病検査の問題の場合を対比させる。

$P_t = P(T)$: ターゲットの事象が真である確率(true)あるいは事前確率。たとえば日本全体での難病の罹患率など。
$P_c = P(T|C)$: 判定が真の条件下で事象が真である確率(under condition)あるいは事後確率。たとえば検査結果が陽性の場合に難病に罹患している確率。
$P_h = P(C|T)$: 事象が真の場合に判定が真となる確率(hit)。たとえば難病に罹患している場合に検査結果が陽性となる確率。
$P_f = P(C| \overline{T})$: 事象が偽なのに判定が正となる確率(fail)。たとえば難病に罹患していないのに検査結果が陽性となる確率。

$P(\overline{T}) = 1 - P(T)$ であることに留意して、これらの記号でベイズの定理を表すと以下の通り。

(3) $\begin{equation*} P_c = \frac{P_h P_t}{P_h P_t + P_F (1 - P_t)} = \frac{1}{1 + \dfrac{P_f}{P_h} \left( \dfrac{1}{P_t} - 1 \right) } \end{equation*}$

いくつかのケース

誤判定を絶対しない場合

(8)において $P_f = 0$ とおくと、 $P_t = 1$ となる。

難病の問題で言えば、「罹患していないときには必ず陰性判定となる(陽性とならない)」という検査に相当し、このような検査の結果が陽性判定なら必ず罹患していることになる。この時、罹患時に陽性となる確率 $P_h$ がいくらかは関係ない。

ただしこの逆は必ずしも保証されておらず、陰性判定でも罹患している場合がある。

完璧な検査

難病の問題で、上記に加えて「罹患しているときは必ず陽性判定となる( $P_h = 1$ 」という条件を加える。

まず、(2)をCの余事象について書くと以下の様になる。

(4) $\begin{eqnarray*} P(T|\overline{C}) &=& \frac{P(\overline{C}|T)P(T)}{P(\overline{C}|T)P(T) + P(\overline{C}|\overline{T})P(\overline{T})} \\ &=& \frac{(1 - P_h) P_t}{(1 - P_h) P_t + (1 - P_f)(1 - P_t)} \\ \end{eqnarray*}$

ここで、 $P(\overline{T} | \overline{C}) = 1 - P(T| \overline{C})$ を考える。これは、陰性反応の時に難病にかかっていない確率に相当する。

(5) $\begin{eqnarray*} P(\overline{T} | \overline{C}) = 1 - P(T| \overline{C}) &=& 1 - \frac{(1 - P_h) P_t}{(1 - P_h) P_t + (1 - P_f)(1 - P_t)} \\ &=& \frac{(1 - P_f)(1 - P_t)}{(1 - P_h) P_t + (1 - P_f)(1 - P_t)} \end{eqnarray*}$

(5)において $P_f = 0$ と置くと、

(6) $\begin{eqnarray*} P(\overline{T} | \overline{C}) = 1 - P(T| \overline{C}) = \frac{(1 - P_t)}{(1 - P_h) P_t + (1 - P_t)} \end{eqnarray*}$

ここで $P_h = 1$ の場合は $P(\overline{T} | \overline{C}) = 1$ となる。

この結果を難病の問題で言えば、検査結果が陰性の場合は必ず罹患していない、ということになる。

以上をまとめると、検査の性質が「罹患していなければ必ず陰性、罹患していれば必ず陽性となる」という場合、検査結果が陽性なら必ず罹患しており、陰性なら罹患していないという、至極当然の結果となる。

意味のない検査

(8)において $P_h = P_f$ と置くと、 $P_c = P_t$ となる。

難病の例で言えば、罹患している時に陽性となる確率と、罹患していないときに陽性となる誤判定率が同じ場合、検査の意味がない(一般的な罹患率でしか罹患の程度がわからない)ということになる。

正しい判定の率が高くても、誤判定率が同じ場合は意味をなさない。

意味のある検査

(8)を $P_f / P_h$ について解く。

(7) $\begin{equation*} \frac{P_f}{P_h} = \frac{\dfrac{1}{P_c} - 1}{\dfrac{1}{P_t}- 1} \end{equation*}$

この式で、 $P_f < P_h$ のとき、 $P_c > P_t$ となり、検査が意味を持つことがわかる。

稀な事象

事前確率 $P_t$ がとても小さい場合を考える。たとえば難病の罹患率が1万人や10万人に1人といった場合。

$P_t \approx 0$ とすると(8)は以下のように変形できる。

(8) $\begin{eqnarray*} P_c &=& \frac{1}{1 + \dfrac{P_f}{P_h} \left( \dfrac{1}{P_t} - 1 \right) } = \frac{P_t}{P_t + \dfrac{P_F}{P_h} (1 - P_t)} \\ &=& \frac{\dfrac{P_h}{P_F} P_t}{1 + \left( \dfrac{P_h}{P_f} - 1 \right) P_t} \approx \frac{P_h}{P_F} P_t} \end{eqnarray*}$

すなわち、事前確率が十分に小さい場合は、検査の正誤判定確率の倍率分だけ事後確率が大きくなる。

たとえば難病検査の問題の場合、 $P_h = 0.9 , \; P_f = 0.2$ として、 $P_t$ が小さくなるほど、事後確率は正誤判定確率の倍率0.9/0.2 = 4.5倍に近づく。

(9) $\begin{eqnarray*} P_t = 0.1 & \rightarrow & P_c = \frac{0.9 \times 0.1}{0.9 \times 0.1 + 0.2 \times 0.9} \approx 0.333 \\ P_t = 0.01 & \rightarrow & P_c = \frac{0.9 \times 0.01}{0.9 \times 0.01 + 0.2 \times 0.99} \approx 0.0435 \\ P_t = 0.001 & \rightarrow & P_c = \frac{0.9 \times 0.001}{0.9 \times 0.01 + 0.2 \times 0.999} \approx 0.00448 \\ P_t = 0.0001 & \rightarrow & P_c = \frac{0.9 \times 0.0001}{0.9 \times 0.01 + 0.2 \times 0.9999} \approx 0.000449 \end{eqnarray*}$