指数分布における両側検定の検出力関数とニュートン法を用いた棄却域の数値計算

数理統計学の検定論で出てくる検出力関数(power function)は一様最強力検定や一様最強力検定の理解にあたって重要な概念です。当記事では指数分布の両側検定の検出力関数を示し、$95$%両側検定における棄却域について多次元ニュートン法を用いて数値計算を行いました。
「数理統計学 統計的推論の基礎(共立出版)」の$10$章の「統計的仮説検定の考え方」の例$10.8$を主に参考に、作成を行いました。

・参考:正規分布の母平均の片側検定・両側検定における検出力関数の描画
https://www.hello-statisticians.com/explain-terms-cat/power_function1.html

指数分布の両側検定の検出力関数

指数分布の確率密度関数

$X \sim \mathrm{Ex}(\lambda)$のとき、確率密度関数を$f(x; \lambda)$とおくと$f(x; \lambda)$は下記のように表される。
$$
\large
\begin{align}
f(x; \lambda) = \lambda e^{-\lambda x}, \, \lambda > 0
\end{align}
$$

指数分布に関しては詳しくは下記などで取り扱った。

両側検定の検出力関数

統計的仮説検定に用いるパラメータ空間$\Theta$を定義するとき、$\theta \in \Theta$の棄却域$C$に関する検出力関数を確率密度関数$f(x; \theta)$を用いて下記のように定める。
$$
\large
\begin{align}
1-\beta(C,\theta) = \int_{C} f(x; \theta) dx
\end{align}
$$

仮説検定では$\Theta$を帰無仮説に対応する$\Theta_0$と対立仮説に対応する$\Theta_1$に分割するので、$\Theta = \Theta_0 \cup \Theta_1$かつ$\phi = \Theta_0 \cap \Theta_1$が成立することも抑えておくと良い。検出力関数に関しては下記で詳しく取りまとめた。

指数分布における両側検定

$X \sim \mathrm{Ex}(\lambda)$のとき、下記のような仮説検定を定義する。
$$
\large
\begin{align}
H_0: \, \lambda = 1 \\
H_1: \, \lambda \neq 1
\end{align}
$$

以下、上記に関して有意水準$\alpha=0.05$の不偏検定の構築を行う。具体的には下記のように定める棄却域$C$の$c_1, c_2$の計算を行う。
$$
\large
\begin{align}
C = (0,c_1) \cup (c_2,\infty)
\end{align}
$$

このとき、検出力関数$1-\beta(C,\lambda)$は下記のように表せる。
$$
\large
\begin{align}
1-\beta(C,\mu) &= \int_{C} f(x; \lambda) dx \\
&= \int_{0}^{c_1} \lambda e^{-\lambda x} dx + \int_{c_2}^{\infty} \lambda e^{-\lambda x} dx \\
&= \left[ -e^{-\lambda x} \right]_{0}^{c_1} + \left[ -e^{-\lambda x} \right]_{c_2}^{\infty} \\
&= -(e^{-c_1 \lambda} – 1) – (0 – e^{-c_2 \lambda}) \\
&= 1 – \exp{(-c_1 \lambda)} + \exp{(-c_2 \lambda)}
\end{align}
$$

ここで$\alpha=0.05$の不偏検定では$\alpha=1-\beta(C,1)=0.05$であるので下記が成立する。
$$
\large
\begin{align}
1-\beta(C,1) = 1 – \exp{(-c_1)} + \exp{(-c_2)} = 0.05 \quad (1)
\end{align}
$$

また、有意水準$\alpha=0.05$の不偏検定であれば検出力関数は$H_0: \, \lambda = 1$より$\lambda=1$で極小値を取るので下記が成立する。
$$
\large
\begin{align}
\frac{d(1-\beta(C,\lambda))}{d \lambda} &= c_1 \exp{(-c_1 \lambda)} – c_2 \exp{(-c_2 \lambda)} \\
\frac{d(1-\beta(C,\lambda))}{d \lambda} \Bigr|_{\lambda=1} &= c_1 \exp{(-c_1)} – c_2 \exp{(-c_2)} = 0 \quad (2)
\end{align}
$$

$(1), (2)$の連立方程式より棄却域に対応する$c_1, c_2$を得ることができる。この連立方程式を解析的に解くことは難しいので、次節では多次元ニュートン法を用いて計算を行う。

棄却域の数値計算

多次元ニュートン法

$$
\large
\begin{align}
\left(\begin{array}{c} x_1^{(n+1)} \\ x_2^{(n+1)} \end{array} \right) &= \left(\begin{array}{c} x_1^{(n)} \\ x_2^{(n)} \end{array} \right) – J(x_1^{(n)},x_2^{(n)})^{-1} \left(\begin{array}{c} f_1(x_1^{(n)},x_2^{(n)}) \\ f_2(x_1^{(n)},x_2^{(n)}) \end{array} \right) \\
J(x_1^{(n)},x_2^{(n)}) &= \left(\begin{array}{cc} \displaystyle \frac{\partial f_1}{\partial x_1} & \displaystyle \frac{\partial f_1}{\partial x_2} \\ \displaystyle \frac{\partial f_2}{\partial x_1} & \displaystyle \frac{\partial f_2}{\partial x_2} \end{array} \middle) \right|_{x_1=x_1^{(n)},x_2=x_2^{(n)}}
\end{align}
$$

二次元のニュートン法は上記のように表すことができる。多次元に拡張する際も同様に考えれば良い。導出は下記で取り扱った。

多次元ニュートン法による棄却域の計算

$$
\large
\begin{align}
\left(\begin{array}{c} f_1(c_1,c_2) \\ f_2(c_1,c_2) \end{array} \right) = \left(\begin{array}{c} 0.95 – \exp(-c_1) + \exp(-c_2) \\ c_1\exp(-c_1) – c_2\exp(-c_2) \end{array} \right) = \left(\begin{array}{c} 0 \\ 0 \end{array} \right)
\end{align}
$$

上記に対し、ヤコビ行列$J(x_1,x_2)$は下記のように計算できる。
$$
\large
\begin{align}
J(c_1,c_2) &= \left(\begin{array}{cc} \displaystyle \frac{\partial f_1}{\partial c_1} & \displaystyle \frac{\partial f_1}{\partial c_2} \\ \displaystyle \frac{\partial f_2}{\partial c_1} & \displaystyle \frac{\partial f_2}{\partial c_2} \end{array} \right) \\
&= \left(\begin{array}{cc} \exp{(-c_1)} & -\exp{(-x_2)} \\ (1-c_1)\exp{(-c_1)} & -(1-c_2)\exp{(-c_2)} \end{array} \right)
\end{align}
$$

上記の連立方程式の解は$(3)$式を元に下記のように計算を行うことができる。

import numpy as np

def calc_f1(c_1, c_2):
    return 0.95 - np.exp(-c_1) + np.exp(-c_2)

def calc_f2(c_1, c_2):
    return c_1*np.exp(-c_1) - c_2*np.exp(-c_2)

c = np.array([[0.7],[1.5]])

for i in range(10):
    J = np.array([[np.exp(-c[0,0]), -np.exp(-c[1,0])], [(1.-c[0,0])*np.exp(-c[0,0]), (c[1,0]-1.)*np.exp(-c[1,0])]])
    c = c - np.linalg.solve(J, np.array([[calc_f1(c[0,0],c[1,0])], [calc_f2(c[0,0],c[1,0])]]))

print(c)

・実行結果

print(calc_f1(0.04236333, 4.76516825))
print(calc_f2(0.04236333, 4.76516825))

上記の結果が適切であることは、下記を実行することで確認することができる。

print(calc_f1(0.04236333, 4.76516825))
print(calc_f2(0.04236333, 4.76516825))

・実行結果

-3.30993582195e-09
-3.0646387858e-09