サンプルサイズの決め方5章 １つの母平均の検定―母分散が未知の場合

はじめに

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

１つの母平均の検定―母分散が未知の場合

母分散が未知の場合の一つの母平均の検定をする手順は書籍では以下のように紹介されています。

• 手順1 仮説の設定

  • 帰無仮説 $H₀: μ = μ₀$ ($μ₀$は指定された値)
  • 対立仮説:
    • (1) $H₁: μ ≠ μ₀$ (両側検定)
    • (2) $H₁: μ > μ₀$ (右片側検定)
    • (3) $H₁: μ < μ₀$ (左片側検定)

• 手順2 有意水準の設定通常は α = 0.05 を使用します。

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

• 手順3 棄却域の設定

  • 対立仮説(1)の場合: $R: |t₀| ≥ t(n-1, α)$
  • 対立仮説(2)の場合: $R: t₀ ≥ t(n-1, 2α)$
  • 対立仮説(3)の場合: $R: t₀ ≤ -t(n-1, 2α)$

ここで、$t(n-1, α)$は$t$分布の上側100α%点です。

• 手順4 検定統計量の計算
  • データ$x₁, x₂, ..., x_n$をとり、検定統計量$t₀$の値を計算:
  • $t₀ = \frac{(x̄ - μ₀)}{\sqrt{(V / n)}}$
  • $V = \frac{1}{(n-1)} Σ(xᵢ - x̄)²$
  • $\phi = n-1$
• 手順5 判定
  • 計算された$t₀$が手順3で設定した棄却域$R$にあれば、有意と判定し、帰無仮説$H₀$を棄却します。

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

import numpy as np
from scipy import stats

def test_population_mean(data, mu0, alpha=0.05, test_type='two-sided'):
    """
    パラメータ:
    data: サンプルデータ
    mu0 (float): 帰無仮説の母平均値
    alpha (float): 有意水準（デフォルトは0.05）
    test_type (str): 検定の種類（'two-sided'：両側検定, 'greater'：右片側検定, 'less'：左片側検定）
    """
    # データをnumpy配列に変換
    data = np.array(data)
    
    # サンプル統計量の計算
    n = len(data)
    x_bar = np.mean(data)
    s = np.std(data, ddof=1)  # 標本標準偏差
    
    # t統計量の計算
    t0 = (x_bar - mu0) / (s / np.sqrt(n))
    
    # 自由度
    df = n - 1
    
    # 対立仮説に基づいて臨界値を決定
    if test_type == 'two-sided':
        cv = stats.t.ppf(1 - alpha/2, df)
        reject = abs(t0) >= cv
    elif test_type == 'greater':
        cv = stats.t.ppf(1 - alpha, df)
        reject = t0 >= cv
    elif test_type == 'less':
        cv = -stats.t.ppf(1 - alpha, df)
        reject = t0 <= cv
    else:
        raise ValueError("無効な検定の種類です。'two-sided'、'greater'、または'less'を選択してください。")
    
    # 結果の準備
    results = {
        "t統計量": t0,
        "自由度": df,
        "臨界値": cv,
        "帰無仮説の棄却": reject,
    }
    
    return results

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

$$\begin{aligned} &\text{問題5.1} \\ &n = 10 \text{ の下記のデータに基づいて、} H_0：\mu = \mu_0 \: (\mu_0 = 5.0)\text{、} H_1：\mu \neq \mu_0 \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}$$

# 問題5.1
data_5_1 = [6.2, 4.8, 7.3, 5.5, 6.5, 4.9, 6.8, 7.9, 6.6, 7.3]
result_5_1 = test_population_mean(data_5_1, mu0=5.0, alpha=0.05, test_type='two-sided')

print("問題5.1の結果:")
print(f"t統計量: {result_5_1['t統計量']:.4f}")
print(f"臨界値: {result_5_1['臨界値']:.4f}")
print(f"帰無仮説を棄却するか: {'はい' if result_5_1['帰無仮説の棄却'] else 'いいえ'}")

結果

問題5.1の結果:
t統計量: 4.1958
臨界値: 2.2622
帰無仮説を棄却するか: はい

$t$統計量(4.1958)が臨界値(2.2622)を超えているため、帰無仮説を棄却します。つまり、母平均が5.0と有意に異なると結論づけられます。

$$\begin{aligned} &\text{問題5.2} \\ &n = 10 \text{ の下記のデータに基づいて、} H_0：\mu = \mu_0 \: (\mu_0 = 10.0)\text{、} H_1：\mu \neq \mu_0 \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}$$

# 問題5.2
data_5_2 = [10.8, 11.2, 9.7, 9.9, 12.0, 9.6, 10.5, 10.7, 10.1]
result_5_2 = test_population_mean(data_5_2, mu0=10.0, alpha=0.05, test_type='greater')

print("問題5.2の結果:")
print(f"t統計量: {result_5_2['t統計量']:.4f}")
print(f"臨界値: {result_5_2['臨界値']:.4f}")
print(f"帰無仮説を棄却するか: {'はい' if result_5_2['帰無仮説の棄却'] else 'いいえ'}")

結果

問題5.2の結果:
t統計量: 1.9285
臨界値: 1.8595
帰無仮説を棄却するか: はい

$t$統計量(1.9285)が臨界値(1.8595)をわずかに超えているため、帰無仮説を棄却します。つまり、母平均が10.0より有意に大きいと結論づけられますが、その差は小さいかもしれません。

$$\begin{aligned} &\text{問題5.3} \\ &n = 10 \text{ の下記のデータに基づいて、} H_0：\mu = \mu_0 \: (\mu_0 = 20.0)\text{、} H_1：\mu \neq \mu_0 \text{ を} \\ &\text{有意水準 } 5\% \text{ で検定せよ。} \\ &\text{データ：} 21,19,16,19, 22,18,20,21 \end{aligned}$$

# 問題5.3
data_5_3 = [21, 19, 16, 19, 22, 18, 20, 21]
result_5_3 = test_population_mean(data_5_3, mu0=20.0, alpha=0.05, test_type='less')

print("問題5.3の結果:")
print(f"t統計量: {result_5_3['t統計量']:.4f}")
print(f"臨界値: {result_5_3['臨界値']:.4f}")
print(f"帰無仮説を棄却するか: {'はい' if result_5_3['帰無仮説の棄却'] else 'いいえ'}")

結果

問題5.3の結果:
t統計量: -0.7338
臨界値: -1.8946
帰無仮説を棄却するか: いいえ

$t$統計量(-0.7338)が臨界値(-1.8946)を超えていないため、帰無仮説を棄却できません。つまり、母平均が20.0より小さいとは言えません。

検出力の計算

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

$$\begin{aligned} \text{両側検定 }H_1: \mu \neq \mu_0 \text{}\\[1em] &Pr(t_0 \leq -t(\phi,\alpha)) + Pr(t_0 \geq t(\phi,\alpha))\\ \text{片側検定 }H_1: \mu > \mu_0\text{} \\[1em] &Pr(t_0 \geq t(\phi, 2\alpha)) \\ \text{片側検定 }H_1: \mu < \mu_0\text{} \\ &Pr(t_0 \leq -t(\phi,2\alpha)) \end{aligned}$$

$t_0$は自由度$\phi=n-1$、非心パラメータ$\lambda=\sqrt{n}Δ$の非心$t$分布$t'(\phi, \lambda)$に従う。ただし、

$$Δ = \frac{\mu-\mu_0}{\sigma}$$

となります。しかし、非心$t$分布の分布関数の計算は複雑であるため、以下のような近似式を使います。

$$Pr(t' \leq w) \approx Pr[\mu \leq \frac{w\{1-1/(4\phi)\}-\lambda}{\sqrt{1+w^2/(2\phi)}}]$$

$t'$は非心t分布に従う確率変数、$w$は任意の値、$\phi$は自由度、$\lambda$は非心パラメータです。右辺は正規分布の確率を表しており、非心t分布を正規分布で近似しています。

上記の検出力の計算を行うpythonプログラムは以下のとおりです。

import numpy as np
from scipy import stats

def calculate_power(n, mu0, mu1, sigma, alpha=0.05, test_type='two-sided'):
    """
    パラメータ:
    n (int): サンプルサイズ
    mu0 (float): 帰無仮説の母平均
    mu1 (float): 対立仮説の母平均
    sigma (float): 母標準偏差
    alpha (float): 有意水準（デフォルトは0.05）
    test_type (str): 検定の種類（'two-sided', 'greater', 'less'）
    """
    # 自由度
    phi = n - 1
    
    # 非心パラメータ
    delta = (mu1 - mu0) / sigma
    lambda_param = np.sqrt(n) * delta
    
    # 臨界値の計算
    if test_type == 'two-sided':
        cv = stats.t.ppf(1 - alpha/2, phi)
    elif test_type in ['greater', 'less']:
        cv = stats.t.ppf(1 - alpha, phi)
    else:
        raise ValueError("無効な検定タイプです。'two-sided', 'greater', 'less'のいずれかを指定してください。")

    # 非心t分布の正規分布による近似
    def noncentral_t_approx(w):
        return stats.norm.cdf((w * (1 - 1/(4*phi)) - lambda_param) / np.sqrt(1 + w**2/(2*phi)))

    # 検出力の計算
    if test_type == 'two-sided':
        power = 1 - noncentral_t_approx(cv) + noncentral_t_approx(-cv)
    elif test_type == 'greater':
        power = 1 - noncentral_t_approx(cv)
    else:  # 'less'
        power = noncentral_t_approx(-cv)

    return power

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

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

# 問題 5.4
print("問題 5.4:")
# (1) n=10, Δ=0.5
power_5_4_1 = calculate_power(n=10, mu0=5.0, mu1=5.0 + 0.5, sigma=1, alpha=0.05, test_type='two-sided')
print(f"(1) n=10, Δ=0.5の場合の検出力: {power_5_4_1:.4f}")

# (2) n=20, Δ=0.5
power_5_4_2 = calculate_power(n=20, mu0=5.0, mu1=5.0 + 0.5, sigma=1, alpha=0.05, test_type='two-sided')
print(f"(2) n=20, Δ=0.5の場合の検出力: {power_5_4_2:.4f}")

結果

問題 5.4:
(1) n=10, Δ=0.5の場合の検出力: 0.2931
(2) n=20, Δ=0.5の場合の検出力: 0.5642

サンプルサイズを10から20に増やすことで、検出力が0.2931から0.5642に大幅に上昇しました。

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

print("\n問題 5.5:")
# (1) n=9, Δ=1.0
power_5_5_1 = calculate_power(n=9, mu0=10.0, mu1=10.0 + 1.0, sigma=1, alpha=0.05, test_type='greater')
print(f"(1) n=9, Δ=1.0の場合の検出力: {power_5_5_1:.4f}")

# (2) n=9, Δ=1.5
power_5_5_2 = calculate_power(n=9, mu0=10.0, mu1=10.0 + 1.5, sigma=1, alpha=0.05, test_type='greater')
print(f"(2) n=9, Δ=1.5の場合の検出力: {power_5_5_2:.4f}")

結果

問題 5.5:
(1) n=9, Δ=1.0の場合の検出力: 0.8615
(2) n=9, Δ=1.5の場合の検出力: 0.9928

$Δ$を1.0から1.5に増やすことで、検出力が0.8615から0.9928に上昇しました。

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

print("\n問題 5.6:")
# (1) n=8, Δ=-1.0
power_5_6_1 = calculate_power(n=8, mu0=20.0, mu1=20.0 - 1.0, sigma=1, alpha=0.05, test_type='less')
print(f"(1) n=8, Δ=-1.0の場合の検出力: {power_5_6_1:.4f}")

# (2) n=16, Δ=-1.0
power_5_6_2 = calculate_power(n=16, mu0=20.0, mu1=20.0 - 1.0, sigma=1, alpha=0.05, test_type='less')
print(f"(2) n=16, Δ=-1.0の場合の検出力: {power_5_6_2:.4f}")

結果

問題 5.6:
(1) n=8, Δ=-1.0の場合の検出力: 0.8142
(2) n=16, Δ=-1.0の場合の検出力: 0.9849

サンプルサイズを8から16に増やすことで、検出力が0.8142から0.9849に上昇しました。

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

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

$$\begin{align*} & \text{両側検定 } (H_0: \mu = \mu_0, \quad H_1: \mu \neq \mu_0) \\[0.5em] & \quad n \approx \left(\frac{z_{\alpha/2} - z_{1-\beta}}{\Delta_{0}}\right)^2 + \frac{z^{2}_{\alpha/2}}{2} \\[2em] & \text{片側検定 } (H_0: \mu = \mu_0, \quad H_1: \mu > \mu_0) \\[0.5em] & \quad n \approx \left(\frac{z_{\alpha} - z_{1-\beta}}{\Delta_{0}}\right)^2 + \frac{z^{2}_{\alpha}}{2} \\[2em] & \text{片側検定 } (H_0: \mu = \mu_0, \quad H_1: \mu < \mu_0) \\[0.5em] & \quad n \approx \left(\frac{z_{\alpha} - z_{1-\beta}}{\Delta_{0}}\right)^2 + \frac{z^{2}_{\alpha}}{2} \\[2em] \end{align*}$$

上記の式通りにサンプルサイズを計算する関数をpythonで実装しました。
この関数を用いて練習問題を解いていきます。

import math
import numpy as np
from scipy import stats

def calculate_sample_size(alpha, beta, delta_0, test_type='two-sided'):
    if test_type == 'two-sided':
        z_alpha = stats.norm.ppf(1 - alpha/2)
    else:
        z_alpha = stats.norm.ppf(1 - alpha)
    z_beta = stats.norm.ppf(1 - beta)
    
    n = ((z_alpha - z_beta) / delta_0) ** 2
    
    n += z_alpha ** 2 / 2
    
    return math.ceil(n)

def find_optimal_sample_size(alpha, target_power, delta, test_type='two-sided'):
    beta = 1 - target_power
    n = calculate_sample_size(alpha, beta, delta, test_type)
    
    while calculate_power(n, 0, delta, 1, alpha, test_type) < target_power:
        n += 1
    
    while calculate_power(n-1, 0, delta, 1, alpha, test_type) >= target_power:
        n -= 1
    
    return n

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

# 問題5.7
alpha = 0.05    # 有意水準
mu0 = 5.0       # 帰無仮説の母平均
delta = 0.8     # Δ = σ/σ0
test_type = 'two-sided'

# (1) 検出力 0.90
target_power_1 = 0.90
n_1 = find_optimal_sample_size(alpha, target_power_1, delta, test_type)
actual_power_1 = calculate_power(n_1, mu0, mu0 + delta, 1, alpha, test_type)
print(f"(1) 検出力 0.90 のときの必要なサンプルサイズ: n = {n_1}")
print(f"   実際の検出力: {actual_power_1:.6f}")

# (2) 検出力 0.95
target_power_2 = 0.95
n_2 = find_optimal_sample_size(alpha, target_power_2, delta, test_type)
actual_power_2 = calculate_power(n_2, mu0, mu0 + delta, 1, alpha, test_type)
print(f"(2) 検出力 0.95 のときの必要なサンプルサイズ: n = {n_2}")
print(f"   実際の検出力: {actual_power_2:.6f}")

結果

(1) 検出力 0.90 のときの必要なサンプルサイズ: n = 19
   実際の検出力: 0.909201
(2) 検出力 0.95 のときの必要なサンプルサイズ: n = 23
   実際の検出力: 0.955900

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

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

# 問題5.8
alpha = 0.05    # 有意水準
mu0 = 10.0       # 帰無仮説の母平均
target_power = 0.90 # 検出力
test_type = 'greater'

# (1) 
delta_1 = 1.0
n_1 = find_optimal_sample_size(alpha, target_power, delta_1, test_type)
actual_power_1 = calculate_power(n_1, mu0, mu0 + delta_1, 1, alpha, test_type)
print(f"(1) 検出力 0.90 のときの必要なサンプルサイズ: n = {n_1}")
print(f"   実際の検出力: {actual_power_1:.6f}")

# (2) 検出力 0.90
delta_2 = 0.5
n_2 = find_optimal_sample_size(alpha, target_power, delta_2, test_type)
actual_power_2 = calculate_power(n_2, mu0, mu0 + delta_2, 1, alpha, test_type)
print(f"(2) 検出力 0.95 のときの必要なサンプルサイズ: n = {n_2}")
print(f"   実際の検出力: {actual_power_2:.6f}")

結果

(1) 検出力 0.90 のときの必要なサンプルサイズ: n = 11
   実際の検出力: 0.924502
(2) 検出力 0.95 のときの必要なサンプルサイズ: n = 36
   実際の検出力: 0.902569

$Δ$を1.0から0.5に下げると、同じ検出力を得るために必要なサンプルサイズが11から36に大幅に増加しました。

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

# 問題5.9
alpha = 0.05    # 有意水準
mu0 = 20.0       # 帰無仮説の母平均
test_type = 'less'

# (1) 検出力 0.90、Δ=-1.2
delta_1 = -1.2
target_power_1 = 0.90 # 検出力
n_1 = find_optimal_sample_size(alpha, target_power_1, delta_1, test_type)
actual_power_1 = calculate_power(n_1, mu0, mu0 + delta_1, 1, alpha, test_type)
print(f"(1) 検出力 0.90 のときの必要なサンプルサイズ: n = {n_1}")
print(f"   実際の検出力: {actual_power_1:.6f}")

# (2) 検出力 0.95、Δ=-1.5
delta_2 = -1.5
target_power_2  = 0.95
n_2 = find_optimal_sample_size(alpha, target_power_2, delta_2, test_type)
actual_power_2 = calculate_power(n_2, mu0, mu0 + delta_2, 1, alpha, test_type)
print(f"(2) 検出力 0.95 のときの必要なサンプルサイズ: n = {n_2}")
print(f"   実際の検出力: {actual_power_2:.6f}")

結果

(1) 検出力 0.90 のときの必要なサンプルサイズ: n = 8
   実際の検出力: 0.918970
(2) 検出力 0.95 のときの必要なサンプルサイズ: n = 7
   実際の検出力: 0.966904

$Δ$の絶対値を大きくすると（-1.2から-1.5）、同じ検出力を得るために必要なサンプルサイズが小さくなりました（8から7）。

最後に

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

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

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