ホーム » 最強力検定
「最強力検定」カテゴリーアーカイブ
2015 Q2(5)
前問の検定がネイマン-ピアソンの基本定理に基づき一様最強力検定であることを示しました。
コード
有意水準α=0.05のもと、検出力を理論値とシミュレーションで描画し比較します。
# 2015 Q2(5) 2024.12.11
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
# シミュレーション設定
alpha = 0.05 # 有意水準
z_alpha = norm.ppf(1 - alpha) # 臨界値 z_alpha
num_simulations = 10000 # シミュレーション回数
n = 10 # 標本サイズ
mu_values = np.linspace(0, 2, 10) # 真の平均 μ の範囲
# 臨界値の計算
critical_value = z_alpha / np.sqrt(n) # 棄却域の閾値
# 検出力を計算する関数
def compute_power_for_mu(mu, n, z_alpha):
return 1 - norm.cdf(z_alpha - mu * np.sqrt(n))
# シミュレーションで検出力を計算
empirical_powers = []
for mu in mu_values:
samples = np.random.normal(mu, 1 / np.sqrt(n), num_simulations) # 標本生成
rejection_rate = np.mean(samples > critical_value) # 棄却割合
empirical_powers.append(rejection_rate)
# 理論的な検出力を計算
theoretical_powers = [compute_power_for_mu(mu, n, z_alpha) for mu in mu_values]
# グラフの描画
plt.figure(figsize=(10, 6))
plt.plot(mu_values, theoretical_powers, label="理論的検出力", linestyle="--", color="blue")
plt.scatter(mu_values, empirical_powers, label="シミュレーションによる検出力", color="red", zorder=5)
plt.axhline(y=0.05, color="green", linestyle="--", label="有意水準 $\\alpha = 0.05$")
# グラフ設定
plt.xlabel("真の平均 $\\mu$", fontsize=12)
plt.ylabel("検出力 (Power)", fontsize=12)
plt.title("理論値とシミュレーションによる検出力の比較", fontsize=14)
plt.legend()
plt.grid(True)
plt.show()
検出力は理論値とシミュレーションでよく一致し、検定が正しく機能していることを確認しました。また、真の平均μ>0に対して検出力が単調増加することが分かり、この検定が一様最強力検定の特徴を備えていることが確認できました。
2019 Q4(4)
コーシー分布の検定で、与えられた棄却域が最強力検定になることを示しました。
コード
棄却域をR={x:1<x<3}にする検定が最強力検定になるのか確認するために、尤度比λ(x)のグラフを描画します。
# 2019 Q4(4) 2024.9.26
import numpy as np
import matplotlib.pyplot as plt
# 尤度比 λ(x) の計算
def likelihood_ratio(x):
return (1 + x**2) / (1 + (x - 1)**2)
# x の範囲を設定
x_values = np.linspace(-5, 5, 500)
# λ(x) の値を計算
lambda_values = likelihood_ratio(x_values)
# グラフを描画
plt.figure(figsize=(10, 6))
plt.plot(x_values, lambda_values, label=r'尤度比 $\lambda(x) = \frac{1 + x^2}{1 + (x - 1)^2}$', color='b')
# 棄却域 1 < x < 3 を塗りつぶす（透明度を追加）
plt.axvspan(1, 3, color='yellow', alpha=0.3, label=r'棄却域 $1 < x < 3$')
# λ(x) >= 2 の部分を赤で塗る
plt.fill_between(x_values, lambda_values, 2, where=(lambda_values >= 2), color='red', alpha=0.3, label=r'$\lambda(x) \geq 2$')
# グラフの装飾
plt.axhline(2, color='green', linestyle='--', label=r'閾値 $\lambda(x) = 2$')
plt.axhline(0, color='black',linewidth=0.5)
plt.axvline(0, color='black',linewidth=0.5)
plt.title(r'尤度比 $\lambda(x)$ と棄却域 $1 < x < 3$, $\lambda(x) \geq 2$ の範囲')
plt.xlabel('x')
plt.ylabel(r'$\lambda(x)$')
plt.grid(True)
plt.legend()
# グラフを表示
plt.show()
棄却域は尤度比λ(x)>c (ここではc=2)で定義できるので、これは最強力検定と言えます。
ところでx->-∞で尤度比λ(x)はどこに収束するのでしょうか。
import numpy as np
import matplotlib.pyplot as plt
# 尤度比 λ(x) の計算
def likelihood_ratio(x):
return (1 + x**2) / (1 + (x - 1)**2)
# x の範囲を設定（左側を広げる）
x_values = np.linspace(-50, 5, 500)
# λ(x) の値を計算
lambda_values = likelihood_ratio(x_values)
# グラフを描画
plt.figure(figsize=(10, 6))
plt.plot(x_values, lambda_values, label=r'尤度比 $\lambda(x) = \frac{1 + x^2}{1 + (x - 1)^2}$', color='b')
# 閾値を視覚化
plt.axhline(1, color='purple', linestyle='--', label=r'極限 $\lambda(x) \to 1$ as $x \to -\infty$')
plt.axhline(0, color='black',linewidth=0.5)
plt.axvline(0, color='black',linewidth=0.5)
plt.title(r'尤度比 $\lambda(x)$ の概形 ($x \to -\infty$ の範囲)')
plt.xlabel('x')
plt.ylabel(r'$\lambda(x)$')
plt.grid(True)
plt.legend()
# グラフを表示
plt.show()
1.0に収束するようです。
よって、cは1以上で、且つ尤度比λ(x)が最大となるx = 1.61803398874989(黄金比)が棄却域に含まれていれば最強力検定になります。