サンプルサイズの決め方4章 １つの母分散の検定

はじめに

このブログ記事では、永田靖著『サンプルサイズの決め方』（朝倉書店）の4章を参考に、母分散の検定、検出力の計算方法、適切なサンプルサイズの設計方法について解説します。これら統計学の基本的な概念をPythonを用いて実装しながら確認していきます。

1つの母分散の検定

母分散を検定する手順は書籍では以下のように紹介されています。

• 手順1 仮説の設定

  • 帰無仮説 $H_0: \sigma^2 = \sigma_0^2$ ($\sigma_0^2$は指定された値)
  • 対立仮説:
    • (1)$H_1: \sigma^2 \neq \sigma_0^2$
    • (2)$H_1: \sigma^2 > \sigma_0^2$
    • (3)$H_1: \sigma^2 < \sigma_0^2$

• 手順2 有意水準の設定

  • 有意水準 $\alpha$ の設定 (通常は $\alpha = 0.05$ とする)

• 手順3 棄却域の設定

  • 棄却域 $R$ の設定:
    • 対立仮説 (1)の場合: $R: \chi_0^2 \leq \chi^2(\phi, 1 - \alpha/2)$ または $\chi_0^2 \geq \chi^2(\phi, \alpha/2)$
    • 対立仮説 (2)の場合: $R: \chi_0^2 \geq \chi^2(\phi, \alpha)$
    • 対立仮説 (3)の場合: $R: \chi_0^2 \leq \chi^2(\phi, 1 - \alpha)$

• 手順4 検定統計量の計算

  • データ $x_1, x_2, \ldots, x_n$ をとり、検定統計量 $\chi_0^2$ の値を計算する。
  • $\chi_0^2 = S / \sigma_0^2$
  • $S = \sum_{i=1}^n (x_i - \bar{x})^2 = \sum_{i=1}^n x_i^2 - (\sum_{i=1}^n x_i)^2 / n$
  • $\phi = n - 1$

• 手順5 判定

  • $\chi_0^2$ が手順3で設定した棄却域 $R$ にあれば有意と判定し、$H_0$ を棄却する。

この検定手順に則って検定を行う関数をpythonで実装しました。

def variance_test(data, sigma0_sq, alpha=0.05, test_type='two-sided'):
    """
    data: サンプルデータ
    sigma0_sq: 帰無仮説の分散
    alpha: 有意水準 (デフォルト 0.05)
    test_type: 'two-sided'（両側検定）, 'greater'（右片側検定） or 'less'（左片側検定）
    :return: 検定統計量, 臨界値, 検定結果
    """
    n = len(data)
    df = n - 1
    sample_var = np.var(data, ddof=0)  # 標本分散を計算（ddof=0）
    chi_sq_stat = n * sample_var / sigma0_sq  # 検定統計量

    if test_type == 'two-sided':
        chi_sq_crit_lower = stats.chi2.ppf(alpha/2, df)
        chi_sq_crit_upper = stats.chi2.ppf(1 - alpha/2, df)
        reject = chi_sq_stat <= chi_sq_crit_lower or chi_sq_stat >= chi_sq_crit_upper
        critical_values = (chi_sq_crit_lower, chi_sq_crit_upper)
    elif test_type == 'greater':
        chi_sq_crit = stats.chi2.ppf(1 - alpha, df)
        reject = chi_sq_stat >= chi_sq_crit
        critical_values = chi_sq_crit
    elif test_type == 'less':
        chi_sq_crit = stats.chi2.ppf(alpha, df)
        reject = chi_sq_stat <= chi_sq_crit
        critical_values = chi_sq_crit
    else:
        raise ValueError("test_type must be 'two-sided', 'greater' or 'less'")

    return chi_sq_stat, critical_values, reject

この関数を用いて練習問題を解いていきます。

$$\begin{aligned} &\text{問題4.1} \\ &n = 10 \text{ の下記のデータに基づいて、} H_0：\sigma^2 = \sigma_0^2 \: (\sigma_0^2 = 2.0^2)\text{、} H_1：\sigma^2 \neq \sigma_0^2 \text{ を} \\ &\text{有意水準 } 5\% \text{ で検定せよ。} \\ &\text{データ：} 6.2, 4.8, 7.3, 5.5, 6.5, 4.9, 6.8, 7.9, 6.6, 7.3 \end{aligned}$$

# 問題4.1
data_4_1 = [6.2, 4.8, 7.3, 5.5, 6.5, 4.9, 6.8, 7.9, 6.6, 7.3]
sigma0_sq_4_1 = 2.0**2
chi_sq_stat, critical_values, reject = variance_test(data_4_1, sigma0_sq_4_1, alpha=0.05, test_type='two-sided')

print(f"検定統計量: {chi_sq_stat:.4f}")
print(f"臨界値: {critical_values}")
print(f"帰無仮説を棄却するか: {'はい' if reject else 'いいえ'}")

結果

検定統計量: 2.4340
臨界値: (2.7003894999803584, 19.02276779864163)
帰無仮説を棄却するか: はい

検定統計量が棄却域にあるため、帰無仮説が棄却されました。このデータは母分散が2.0²と異なることを示唆しています。

$$\begin{aligned} &\text{問題4.2} \\ &n = 9 \text{ の下記のデータに基づいて、} H_0：\sigma^2 = \sigma_0^2 \: (\sigma_0^2 = 0.6^2)\text{、} H_1：\sigma^2 > \sigma_0^2 \text{ を} \\ &\text{有意水準 } 5\% \text{ で検定せよ。} \\ &\text{データ：} 10.8, 11.2, 9.7, 9.9, 12.0, 9.6, 10.5, 10.7, 10.1 \end{aligned}$$

# 問題4.2
data_4_2 = [10.8, 11.2, 9.7, 9.9, 12.0, 9.6, 10.5, 10.7, 10.1]
sigma0_sq_4_2 = 0.6**2
chi_sq_stat, critical_values, reject = variance_test(data_4_2, sigma0_sq_4_2, alpha=0.05, test_type='greater')

print(f"検定統計量: {chi_sq_stat:.4f}")
print(f"臨界値: {critical_values:.4f}")
print(f"帰無仮説を棄却するか: {'はい' if reject else 'いいえ'}")

結果

検定統計量: 13.4444
臨界値: 15.5073
帰無仮説を棄却するか: いいえ

検定統計量が棄却域にないため、帰無仮説を棄却できません。このデータからは、母分散が0.6²より大きいとは結論づけられません。

$$\begin{aligned} &\text{問題4.3} \\ &n = 8 \text{ の下記のデータに基づいて、} H_0：\sigma^2 = \sigma_0^2 \: (\sigma_0^2 = 4.0^2)\text{、} H_1：\sigma^2 < \sigma_0^2 \text{ を} \\ &\text{有意水準 } 5\% \text{ で検定せよ。} \\ &\text{データ：} 21, 19, 16, 19, 22, 18, 20, 21 \end{aligned}$$

# 問題4.3
data_4_3 = [21, 19, 16, 19, 22, 18, 20, 21]
sigma0_sq_4_3 = 4.0**2
chi_sq_stat, critical_value, reject = variance_test(data_4_3, sigma0_sq_4_3, alpha=0.05, test_type='less')

print(f"検定統計量: {chi_sq_stat:.4f}")
print(f"臨界値: {critical_value:.4f}")
print(f"帰無仮説を棄却するか: {'はい' if reject else 'いいえ'}")

結果

検定統計量: 1.6250
臨界値: 2.1673
帰無仮説を棄却するか: はい

検定統計量が棄却域にあるため、帰無仮説が棄却されました。このデータは母分散が4.0²未満であることを示唆しています。

検出力の計算

検出力を計算する手順は書籍では以下のように紹介されています。

$$\begin{aligned} \text{両側検定 }H_1: \sigma^2 \neq \sigma_0^2 \text{}\\[1em] &Pr\left(\chi^2 \leq \frac{\chi^2(n-1,1-\alpha/2)}{\Delta^2}\right) + Pr\left(\chi^2 \geq \frac{\chi^2(n-1,\alpha/2)}{\Delta^2}\right)\\ \text{片側検定 }H_1: \sigma^2 > \sigma_0^2\text{} \\[1em] &Pr\left(\chi^2 \geq \frac{\chi^2(n-1,\alpha)}{\Delta^2}\right) \\ \text{片側検定 }H_1: \sigma^2 < \sigma_0^2\text{} \\ &Pr\left(\chi^2 \leq \frac{\chi^2(n-1,1-\alpha)}{\Delta^2}\right) \end{aligned}$$

以下のプログラムは検出力の計算をpythonで実装したものとなっています。
このcalculate_power関数を用いて、練習問題を解いてきます。

def calculate_power(n, alpha, delta, test_type='two-sided'):
    """
    検出力を計算する関数
    
    n: サンプルサイズ
    alpha: 有意水準
    delta: Δ 
    test_type: 'two-sided', 'greater' or 'less'
    :return: 検出力 (1 - β)
    """
    df = n - 1
    
    if test_type == 'two-sided':
        chi2_lower = stats.chi2.ppf(alpha/2, df)
        chi2_upper = stats.chi2.ppf(1 - alpha/2, df)
        power = stats.chi2.cdf(chi2_lower / delta**2, df) + \
                (1 - stats.chi2.cdf(chi2_upper / delta**2, df))
    elif test_type == 'greater':
        chi2_crit = stats.chi2.ppf(1 - alpha, df)
        power = 1 - stats.chi2.cdf(chi2_crit / delta**2, df)
    elif test_type == 'less':
        chi2_crit = stats.chi2.ppf(alpha, df)
        power = stats.chi2.cdf(chi2_crit / delta**2, df)
    else:
        raise ValueError("test_type must be 'two-sided', 'greater' or 'less'")
    
    return power

$$\begin{aligned} &\text{問題4.4} \\ &\text{帰無仮説と対立仮説を } H_0：\sigma^2 = \sigma^2_0 \: (\sigma^2_0 = 2.0^2)\text{、} H_1：\sigma^2 ≠ \sigma^2_0 \text{ と設定し、有意水準 } 5\% \text{ で検定する。} \\ &\text{(1) } n = 10\text{、}\Delta = \sigma/\sigma_0 = 0.5 \text{ の場合に検出力を求めよ。} \\ &\text{(2) } n = 20\text{、}\Delta = \sigma/\sigma_0 = 0.5 \text{ の場合に検出力を求めよ。} \end{aligned}$$

print("問題4.4の結果:")
alpha_4_4 = 0.05
delta_4_4 = 0.5

# (1) n = 10
n_4_4_1 = 10
power_4_4_1 = calculate_power(n_4_4_1, alpha_4_4, delta_4_4, 'two-sided')
print(f"(1) n = 10 の場合の検出力: {power_4_4_1:.4f}")

# (2) n = 20
n_4_4_2 = 20
power_4_4_2 = calculate_power(n_4_4_2, alpha_4_4, delta_4_4, 'two-sided')
print(f"(2) n = 20 の場合の検出力: {power_4_4_2:.4f}")

結果

問題4.4の結果:
(1) n = 10 の場合の検出力: 0.7104
(2) n = 20 の場合の検出力: 0.9883

サンプルサイズを10から20に増やすことで、検出力が0.7104から0.9883へと向上しています。これは、より大きなサンプルサイズが検出力が高くなることを示しています。

$$\begin{aligned} &\text{問題4.5} \\ &\text{帰無仮説と対立仮説を } H_0：\sigma^2 = \sigma^2_0 \: (\sigma^2_0 = 0.6^2)\text{、} H_1：\sigma^2 > \sigma^2_0 \text{ と設定し、有意水準 } 5\% \text{ で検定する。} \\ &\text{(1) } n = 9\text{、}\Delta = \sigma/\sigma_0 = 1.5 \text{ の場合に検出力を求めよ。} \\ &\text{(2) } n = 9\text{、}\Delta = \sigma/\sigma_0 = 2.0 \text{ の場合に検出力を求めよ。} \end{aligned}$$

print("問題4.5の結果:")
alpha_4_5 = 0.05
n = 9

# (1) delta = 1.5
delta_4_5_1 = 1.5
power_4_4_1 = calculate_power(n, alpha_4_5, delta_4_5_1, 'greater')
print(f"(1) delta = 1.5 の場合の検出力: {power_4_4_1:.4f}")

# (2) delta = 2.0
delta_4_5_2 = 2.0
power_4_4_2 = calculate_power(n, alpha_4_5, delta_4_5_2, 'greater')
print(f"(2) delta = 2.0 の場合の検出力: {power_4_4_2:.4f}")

結果

問題4.5の結果:
(1) delta = 1.5 の場合の検出力: 0.5483
(2) delta = 2.0 の場合の検出力: 0.8681

Δを1.5から2.0に増やすことで、検出力が0.5483から0.8681に上昇しています。

$$\begin{aligned} &\text{問題4.6} \\ &\text{帰無仮説と対立仮説を } H_0：\sigma^2 = \sigma^2_0 \: (\sigma^2_0 = 4.0^2)\text{、} H_1：\sigma^2 < \sigma^2_0 \text{ と設定し、有意水準 } 5\% \text{ で検定する。} \\ &\text{(1) } n = 8\text{、}\Delta = \sigma/\sigma_0 = 0.8 \text{ の場合に検出力を求めよ。} \\ &\text{(2) } n = 16\text{、}\Delta = \sigma/\sigma_0 = 0.8 \text{ の場合に検出力を求めよ。} \end{aligned}$$

print("問題4.6の結果:")
alpha_4_6 = 0.05
delta_4_6 = 0.8

# (1) n = 8
n_4_6_1 = 8
power_4_6_1 = calculate_power(n_4_6_1, alpha_4_6, delta_4_6, 'less')
print(f"(1) n = 8 の場合の検出力: {power_4_6_1:.4f}")

# (2) n = 16
n_4_6_2 = 16
power_4_6_2 = calculate_power(n_4_6_2, alpha_4_6, delta_4_6, 'less')
print(f"(2) n = 16 の場合の検出力: {power_4_6_2:.4f}")

結果

問題4.6の結果:
(1) n = 8 の場合の検出力: 0.1529
(2) n = 16 の場合の検出力: 0.2722

サンプルサイズを8から16に倍増させると、検出力は0.1529から0.2722へと増加しています。

サンプルサイズの設計方法

サンプルサイズの設計方法は書籍では以下のように紹介されています。

$$\begin{align*} & \text{両側検定 } (H_0: \sigma^2 = \sigma_0^2, \quad H_1: \sigma^2 \neq \sigma_0^2) \\[0.5em] & \quad \text{(a) } \Delta_{01} > 1 \text{ の場合:} \quad n \approx \frac{1}{2}\left(\frac{z_{\alpha/2} - \Delta_{01}z_{1-\beta}}{\Delta_{01} - 1}\right)^2 + \frac{3}{2} \\[1em] & \quad \text{(b) } \Delta_{02} < 1 \text{ の場合:} \quad n \approx \frac{1}{2}\left(\frac{z_{\alpha/2} - \Delta_{02}z_{1-\beta}}{\Delta_{02} - 1}\right)^2 + \frac{3}{2} \\[2em] & \text{片側検定 } (H_0: \sigma^2 = \sigma_0^2, \quad H_1: \sigma^2 > \sigma_0^2) \\[0.5em] & \quad n \approx \frac{1}{2}\left(\frac{z_{\alpha} - \Delta_0z_{1-\beta}}{\Delta_0 - 1}\right)2 + \frac{3}{2} \\[2em] & \text{片側検定 } (H_0: \sigma^2 = \sigma_0^2, \quad H_1: \sigma^2 < \sigma_0^2) \\[0.5em] & \quad n \approx \frac{1}{2}\left(\frac{z_{\alpha} - \Delta_0z_{1-\beta}}{\Delta_0 - 1}\right)2 + \frac{3}{2} \end{align*}$$

上記の式通りにサンプルサイズを計算する関数をpythonで実装しました。
これと検出力の計算の部分で実装したcalculate_powerを利用して、練習問題を解いていきます。

def calculate_sample_size(alpha, power, delta, test_type='two-sided'):
    z_alpha = stats.norm.ppf(1 - alpha/2) if test_type == 'two-sided' else stats.norm.ppf(1 - alpha)
    z_beta = stats.norm.ppf(power)
    
    if test_type == 'two-sided':
        n1 = 0.5 * ((z_alpha - delta * z_beta) / (delta - 1))**2 + 1.5
        n2 = 0.5 * ((z_alpha - (1/delta) * z_beta) / (1/delta - 1))**2 + 1.5
        return n1, n2
    else:
        n = 0.5 * ((z_alpha - delta * z_beta) / (delta - 1))**2 + 1.5
        return n

def find_optimal_sample_size(alpha, target_power, delta, test_type='two-sided'):
    if test_type == 'two-sided':
        n1, n2 = calculate_sample_size(alpha, 1 - target_power, delta, test_type)
        n1 = max(2, int(np.ceil(n1)))
        n2 = max(2, int(np.ceil(n2)))
        
        while calculate_power(n1, alpha, delta, test_type) > target_power:
            n1 -= 1
        n1 += 1
        
        while calculate_power(n2, alpha, 1/delta, test_type) > target_power:
            n2 -= 1
        n2 += 1
        
        return max(2, n1), max(2, n2)
    else:
        n = max(2, int(np.ceil(calculate_sample_size(alpha, 1 - target_power, delta, test_type))))
        while calculate_power(n, alpha, delta, test_type) > target_power:
            n -= 1
        n += 1
        return max(2, n)

$$\begin{aligned} &\text{問題4.7} \\ &\text{帰無仮説と対立仮説を } H_0：\sigma^2 = \sigma^2_0 \: (\sigma^2_0 = 2.0^2)\text{、} H_1：\sigma^2 \neq \sigma^2_0 \text{ と設定し、}\\ &\text{有意水準 } 5\% \text{ で検定する。 } \\[0.5em] &\text{(1) } \Delta = \sigma/\sigma_0 \geq -1.8 \text{または} \Delta = \sigma/\sigma_0 \leq 5/9\text{ のとき検出力が } 0.90 \text{ となるように} \\ &\quad\text{サンプルサイズ } n \text{ を定めよ。} \\[0.5em] &\text{(2) } \Delta = \sigma/\sigma_0 \geq -1.8 \text{または} \Delta = \sigma/\sigma_0 \leq 5/9\text{ のとき検出力が } 0.95 \text{ となるように} \\ &\quad\text{サンプルサイズ } n \text{ を定めよ。} \\[0.5em] \end{aligned}$$

# 問題4.7
print("問題4.7:")
alpha = 0.05
delta1 = 1.8
delta2 = 5/9

# (1) 検出力 0.90
target_power1 = 0.90
n1_1, n2_1 = find_optimal_sample_size(alpha, target_power1, delta1, 'two-sided')
n1_2, n2_2 = find_optimal_sample_size(alpha, target_power1, delta2, 'two-sided')

print(f"(1) 検出力 {target_power1}:")
print(f"  Δ = {delta1} の場合: n = {n1_1}")
print(f"  Δ = {1/delta1} の場合: n = {n2_1}")

# (2) 検出力 0.95
target_power2 = 0.95
n1_1, n2_1 = find_optimal_sample_size(alpha, target_power2, delta1, 'two-sided')
n1_2, n2_2 = find_optimal_sample_size(alpha, target_power2, delta2, 'two-sided')

print(f"\n(2) 検出力 {target_power2}:")
print(f"  Δ = {delta1} の場合: n = {n1_1}")
print(f"  Δ = {1/delta1} の場合: n = {n2_1}")

結果

問題4.7:
(1) 検出力 0.9:
  Δ = 1.8 の場合: n = 16
  Δ = 0.5555555555555556 の場合: n = 19

(2) 検出力 0.95:
  Δ = 1.8 の場合: n = 21
  Δ = 0.5555555555555556 の場合: n = 22

検出力を0.90から0.95に上げるには、必要なサンプルサイズを増やす必要があることがわかります。

$$\begin{aligned} &\text{問題4.8} \\ &\text{帰無仮説と対立仮説を } H_0：\sigma^2 = \sigma^2_0 \: (\sigma^2_0 = 0.6^2)\text{、} H_1：\sigma^2 > \sigma^2_0 \text{と設定し、}\\ &\text{有意水準 } 5\% \text{ で検定する。 } \\[0.5em] &\text{(1) } \Delta = \sigma/\sigma_0 \geq 2.0 \text{ のとき検出力が } 0.90 \text{ となるようにサンプルサイズ } n \text{ を定めよ。} \\[0.5em] &\text{(2) } \Delta = \sigma/\sigma_0 \geq 2.0 \text{ のとき検出力が } 0.95 \text{ となるようにサンプルサイズ } n \text{ を定めよ。} \\[0.5em] \end{aligned}$$

# 問題4.8
print("\n問題4.8:")
alpha = 0.05
delta = 2.0

# (1) 検出力 0.90
target_power1 = 0.90
n1 = find_optimal_sample_size(alpha, target_power1, delta, 'greater')
print(f"(1) 検出力 {target_power1}:")
print(f"  必要サンプルサイズ: {n1}")

# (2) 検出力 0.95
target_power2 = 0.95
n2 = find_optimal_sample_size(alpha, target_power2, delta, 'greater')
print(f"\n(2) 検出力 {target_power2}:")
print(f"  必要サンプルサイズ: {n2}")

結果

問題4.8:
(1) 検出力 0.9:
  必要サンプルサイズ: 11

(2) 検出力 0.95:
  必要サンプルサイズ: 14

$$\begin{aligned} &\text{問題4.9} \\ &\text{帰無仮説と対立仮説を } H_0：\sigma^2 = \sigma^2_0 \: (\sigma^2_0 = 4.0^2)\text{、} H_1：\sigma^2 < \sigma^2_0 \text{と設定し、}\\ &\text{有意水準 } 5\% \text{ で検定する。 } \\[0.5em] &\text{(1) } \Delta = \sigma/\sigma_0 \leq 1/2 \text{ のとき検出力が } 0.90 \text{ となるようにサンプルサイズ } n \text{ を定めよ。} \\[0.5em] &\text{(2) } \Delta = \sigma/\sigma_0 \leq 1/2 \text{ のとき検出力が } 0.95 \text{ となるようにサンプルサイズ } n \text{ を定めよ。} \\[0.5em] \end{aligned}$$

# 問題4.9
print("\n問題4.9:")
alpha = 0.05
delta = 1/2

# (1) 検出力 0.90
target_power1 = 0.90
n1 = find_optimal_sample_size(alpha, target_power1, delta, 'less')
print(f"(1) 検出力 {target_power1}:")
print(f"  必要サンプルサイズ: {n1}")

# (2) 検出力 0.95
target_power2 = 0.95
n2 = find_optimal_sample_size(alpha, target_power2, delta, 'less')
print(f"\n(2) 検出力 {target_power2}:")
print(f"  必要サンプルサイズ: {n2}")

結果

問題4.9:
(1) 検出力 0.9:
  必要サンプルサイズ: 12

(2) 検出力 0.95:
  必要サンプルサイズ: 14

最後に

本ブログでは、永田靖著『サンプルサイズの決め方』の第4章「１つの母分散の検定」の内容を、Pythonを用いて実践的に解説しました。ここで学んだ主要なポイントは以下の通りです：

1. 母分散の検定の基本概念と手順
2. 検出力の計算方法とその重要性
3. 適切なサンプルサイズの設計方法

次回のブログでは、『サンプルサイズの決め方』の第5章の内容を取り上げる予定です。このブログを通して少しでもサンプルサイズの設計方法について理解が進めば幸いです。