統計学では、推定量 $\hat{\theta}$ の漸近分布（$\hat{\theta}$ が近似的に従う正規分布）がわかっている場合に、その関数 $g(\hat{\theta})$ の漸近分布を知りたいことが頻繁にあります。例えば、確率の推定値 $\hat{p}$ からオッズ $\hat{p}/(1-\hat{p})$ や対数オッズ $\log(\hat{p}/(1-\hat{p}))$ の分散や信頼区間を求めたい場合です。
デルタ法（delta method） は、テイラー展開を用いて $g(\hat{\theta})$ の漸近分布を導出する手法です。簡単な微分計算だけで分散の近似公式が得られるため、統計学のほぼすべての分野で日常的に使用されています。
本記事の内容
- デルタ法の動機と直感
- 1変数デルタ法の導出（テイラー展開 + スラツキーの定理）
- 多変数デルタ法の導出
- 2次デルタ法（バイアス補正項）
- 具体例（対数変換、比率の分散、オッズ比の信頼区間）
- ロジスティック回帰の係数の信頼区間への応用
- Pythonでのデルタ法とブートストラップ法の比較
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
デルタ法の動機
問題設定
$X_1, X_2, \dots, X_n$ が iid で $E[X_i] = \mu$、$\text{Var}(X_i) = \sigma^2$ のとき、中心極限定理により
$$ \sqrt{n}(\bar{X}_n – \mu) \xrightarrow{d} N(0, \sigma^2) $$
あるいは同等に、
$$ \bar{X}_n \overset{\text{approx}}{\sim} N\left(\mu, \frac{\sigma^2}{n}\right) $$
ここで、関数 $g$ を適用した統計量 $g(\bar{X}_n)$ の漸近分布を知りたいとします。例えば、
- $g(x) = x^2$（分散の推定に関連）
- $g(x) = \log x$（対数変換）
- $g(x) = 1/x$（逆数変換）
- $g(x) = e^x$（指数変換）
$g$ が非線形関数の場合、$g(\bar{X}_n)$ の分布は一般に正規分布ではありません。しかし、$n$ が大きいとき $\bar{X}_n$ は $\mu$ の近くに集中するため、$g$ を $\mu$ の周りで線形近似できます。これがデルタ法の核心的なアイデアです。
1変数デルタ法の導出
定理の主張
定理（1変数デルタ法）： 確率変数の列 $\{T_n\}$ が
$$ \sqrt{n}(T_n – \theta) \xrightarrow{d} N(0, \sigma^2) $$
を満たすとする。$g$ が $\theta$ で微分可能で $g'(\theta) \neq 0$ ならば、
$$ \begin{equation} \sqrt{n}(g(T_n) – g(\theta)) \xrightarrow{d} N(0, \sigma^2 [g'(\theta)]^2) \end{equation} $$
あるいは等価に、$g(T_n)$ の漸近分散は
$$ \begin{equation} \text{Var}(g(T_n)) \approx \frac{\sigma^2 [g'(\theta)]^2}{n} \end{equation} $$
テイラー展開を用いた導出
$g(T_n)$ を $\theta$ の周りで1次のテイラー展開します。
$$ g(T_n) = g(\theta) + g'(\theta)(T_n – \theta) + \frac{1}{2}g”(\theta^*)(T_n – \theta)^2 $$
ここで $\theta^*$ は $\theta$ と $T_n$ の間の値です（テイラーの定理の剰余項）。
両辺から $g(\theta)$ を引き、$\sqrt{n}$ を掛けると、
$$ \sqrt{n}(g(T_n) – g(\theta)) = g'(\theta) \cdot \sqrt{n}(T_n – \theta) + \frac{1}{2}g”(\theta^*)\sqrt{n}(T_n – \theta)^2 $$
第2項について評価します。$\sqrt{n}(T_n – \theta) \xrightarrow{d} N(0, \sigma^2)$ より $\sqrt{n}(T_n – \theta) = O_p(1)$ です。また $T_n – \theta = O_p(1/\sqrt{n})$ なので、
$$ \sqrt{n}(T_n – \theta)^2 = \sqrt{n} \cdot O_p(1/n) = O_p(1/\sqrt{n}) \to 0 $$
よって剰余項は確率的にゼロに収束します。
$$ \sqrt{n}(g(T_n) – g(\theta)) = g'(\theta) \cdot \sqrt{n}(T_n – \theta) + o_p(1) $$
スラツキーの定理 により、$o_p(1)$ は漸近分布に影響しないので、
$$ \sqrt{n}(g(T_n) – g(\theta)) \xrightarrow{d} g'(\theta) \cdot N(0, \sigma^2) = N(0, \sigma^2[g'(\theta)]^2) $$
$\blacksquare$
スラツキーの定理の復習
スラツキーの定理： $X_n \xrightarrow{d} X$ かつ $Y_n \xrightarrow{p} c$（定数）ならば、
- $X_n + Y_n \xrightarrow{d} X + c$
- $X_n \cdot Y_n \xrightarrow{d} cX$
上の導出では、$X_n = g'(\theta)\sqrt{n}(T_n – \theta)$ と $Y_n = o_p(1) \xrightarrow{p} 0$ に対して適用しました。
多変数デルタ法の導出
定理の主張
定理（多変数デルタ法）： $p$ 次元の確率変数ベクトルの列 $\{\bm{T}_n\}$ が
$$ \sqrt{n}(\bm{T}_n – \bm{\theta}) \xrightarrow{d} N(\bm{0}, \bm{\Sigma}) $$
を満たすとする。$g: \mathbb{R}^p \to \mathbb{R}$ が $\bm{\theta}$ で微分可能で、勾配ベクトル $\nabla g(\bm{\theta}) \neq \bm{0}$ ならば、
$$ \begin{equation} \sqrt{n}(g(\bm{T}_n) – g(\bm{\theta})) \xrightarrow{d} N\left(0, \nabla g(\bm{\theta})^T \bm{\Sigma}\, \nabla g(\bm{\theta})\right) \end{equation} $$
ここで $\nabla g(\bm{\theta}) = \left(\frac{\partial g}{\partial \theta_1}, \dots, \frac{\partial g}{\partial \theta_p}\right)^T$ は $\bm{\theta}$ での勾配ベクトルです。
導出
多変数テイラー展開の1次近似は、
$$ g(\bm{T}_n) \approx g(\bm{\theta}) + \nabla g(\bm{\theta})^T (\bm{T}_n – \bm{\theta}) $$
両辺から $g(\bm{\theta})$ を引いて $\sqrt{n}$ を掛けると、
$$ \sqrt{n}(g(\bm{T}_n) – g(\bm{\theta})) \approx \nabla g(\bm{\theta})^T \sqrt{n}(\bm{T}_n – \bm{\theta}) $$
$\sqrt{n}(\bm{T}_n – \bm{\theta}) \xrightarrow{d} N(\bm{0}, \bm{\Sigma})$ と、正規分布の線形変換の性質（$\bm{a}^T \bm{Z}$ の分散は $\bm{a}^T\bm{\Sigma}\bm{a}$）より、
$$ \sqrt{n}(g(\bm{T}_n) – g(\bm{\theta})) \xrightarrow{d} N(0, \nabla g(\bm{\theta})^T \bm{\Sigma}\, \nabla g(\bm{\theta})) $$
$\blacksquare$
ベクトル値関数への拡張
$\bm{g}: \mathbb{R}^p \to \mathbb{R}^q$ のとき、ヤコビ行列 $\bm{J} = \frac{\partial \bm{g}}{\partial \bm{\theta}^T}$（$q \times p$ 行列）を用いて、
$$ \sqrt{n}(\bm{g}(\bm{T}_n) – \bm{g}(\bm{\theta})) \xrightarrow{d} N(\bm{0}, \bm{J}\bm{\Sigma}\bm{J}^T) $$
2次デルタ法
バイアスの問題
1次デルタ法は漸近分散の近似としては優れていますが、$g(T_n)$ の期待値について見ると、
$$ E[g(T_n)] \approx g(\theta) $$
しか得られず、バイアスの情報が含まれていません。2次のテイラー展開まで含めると、バイアスの補正が可能になります。
2次デルタ法の導出
$g(T_n)$ を2次まで展開します。
$$ g(T_n) = g(\theta) + g'(\theta)(T_n – \theta) + \frac{1}{2}g”(\theta)(T_n – \theta)^2 + O_p(n^{-3/2}) $$
期待値を取ると、$E[T_n – \theta] \approx 0$（不偏性を仮定）のとき、
$$ \begin{align} E[g(T_n)] &\approx g(\theta) + g'(\theta) \cdot 0 + \frac{1}{2}g”(\theta) \cdot E[(T_n – \theta)^2] \\ &= g(\theta) + \frac{1}{2}g”(\theta) \cdot \text{Var}(T_n) \\ &= g(\theta) + \frac{g”(\theta)\sigma^2}{2n} \end{align} $$
したがって、$g(T_n)$ の バイアス は近似的に
$$ \begin{equation} \text{Bias}[g(T_n)] \approx \frac{g”(\theta)\sigma^2}{2n} \end{equation} $$
2次デルタ法による分散
分散についても2次まで含めると、
$$ \begin{align} \text{Var}(g(T_n)) &\approx [g'(\theta)]^2 \text{Var}(T_n) + \frac{1}{2}[g”(\theta)]^2 \text{Var}((T_n – \theta)^2) \\ &\quad + g'(\theta)g”(\theta)\text{Cov}(T_n – \theta, (T_n – \theta)^2) \end{align} $$
$T_n$ が漸近正規のとき、$(T_n – \theta)^2$ の分散は $O(1/n^2)$ の項なので、主要項は
$$ \text{Var}(g(T_n)) \approx \frac{[g'(\theta)]^2\sigma^2}{n} $$
であり、1次デルタ法と一致します。2次項は $O(1/n^2)$ で、標本サイズが大きいときには影響が小さくなります。
具体例
例1：対数変換
$X_1, \dots, X_n$ が iid で $E[X_i] = \mu > 0$、$\text{Var}(X_i) = \sigma^2$ のとき、$g(x) = \log x$ の場合。
$g'(x) = 1/x$ なので、
$$ \text{Var}(\log \bar{X}_n) \approx \frac{\sigma^2}{n\mu^2} $$
2次デルタ法によるバイアス：$g”(x) = -1/x^2$ なので、
$$ \text{Bias}[\log \bar{X}_n] \approx \frac{-\sigma^2}{2n\mu^2} $$
これは $\log \bar{X}_n$ が $\log \mu$ を過小推定する傾向があることを示しています（ジェンセンの不等式 $E[\log X] \leq \log E[X]$ と整合的です）。
例2：比率の分散
2つの独立な推定量 $\hat{\mu}_1 \sim N(\mu_1, \sigma_1^2/n)$ と $\hat{\mu}_2 \sim N(\mu_2, \sigma_2^2/n)$ の比 $R = \hat{\mu}_1 / \hat{\mu}_2$ の分散を求めます。
$g(\mu_1, \mu_2) = \mu_1 / \mu_2$ のとき、
$$ \frac{\partial g}{\partial \mu_1} = \frac{1}{\mu_2}, \quad \frac{\partial g}{\partial \mu_2} = -\frac{\mu_1}{\mu_2^2} $$
多変数デルタ法により（$\hat{\mu}_1$ と $\hat{\mu}_2$ が独立なので $\bm{\Sigma}$ は対角行列）、
$$ \begin{align} \text{Var}(R) &\approx \left(\frac{1}{\mu_2}\right)^2 \frac{\sigma_1^2}{n} + \left(\frac{\mu_1}{\mu_2^2}\right)^2 \frac{\sigma_2^2}{n} \\ &= \frac{1}{n}\left(\frac{\sigma_1^2}{\mu_2^2} + \frac{\mu_1^2 \sigma_2^2}{\mu_2^4}\right) \\ &= \frac{1}{n\mu_2^2}\left(\sigma_1^2 + R^2\sigma_2^2\right) \end{align} $$
ここで $R = \mu_1/\mu_2$ です。相対誤差の形に整理すると、
$$ \frac{\text{Var}(R)}{R^2} \approx \frac{1}{n}\left(\frac{\sigma_1^2}{\mu_1^2} + \frac{\sigma_2^2}{\mu_2^2}\right) $$
比の相対分散は、分子と分母の変動係数（coefficient of variation）の2乗の和に比例します。
例3：オッズ比の対数の信頼区間
2つのベルヌーイ集団で $\hat{p}_1 \sim N(p_1, p_1(1-p_1)/n_1)$、$\hat{p}_2 \sim N(p_2, p_2(1-p_2)/n_2)$ とします。オッズ比は
$$ \text{OR} = \frac{p_1/(1-p_1)}{p_2/(1-p_2)} $$
対数オッズ比を $g(p_1, p_2) = \log\frac{p_1(1-p_2)}{p_2(1-p_1)}$ とします。
偏微分を計算します。
$$ \frac{\partial g}{\partial p_1} = \frac{1}{p_1} + \frac{1}{1-p_1} = \frac{1}{p_1(1-p_1)} $$
$$ \frac{\partial g}{\partial p_2} = -\frac{1}{p_2} – \frac{1}{1-p_2} = -\frac{1}{p_2(1-p_2)} $$
$\hat{p}_1$ と $\hat{p}_2$ は独立なので、
$$ \begin{align} \text{Var}(\log \widehat{\text{OR}}) &\approx \frac{1}{[p_1(1-p_1)]^2} \cdot \frac{p_1(1-p_1)}{n_1} + \frac{1}{[p_2(1-p_2)]^2} \cdot \frac{p_2(1-p_2)}{n_2} \\ &= \frac{1}{n_1 p_1(1-p_1)} + \frac{1}{n_2 p_2(1-p_2)} \end{align} $$
$2 \times 2$ 分割表で $a, b, c, d$ をセルの度数とすると、$\hat{p}_1 = a/(a+b)$、$n_1 = a+b$ 等を代入して、近似的に
$$ \text{Var}(\log \widehat{\text{OR}}) \approx \frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d} $$
という有名な公式が得られます。対数オッズ比の95%信頼区間は
$$ \log \widehat{\text{OR}} \pm 1.96\sqrt{\frac{1}{a} + \frac{1}{b} + \frac{1}{c} + \frac{1}{d}} $$
指数を取って元のスケールに戻すとオッズ比の信頼区間が得られます。
例4：ロジスティック回帰の係数の信頼区間
ロジスティック回帰 $\log\frac{P(Y=1|\bm{x})}{1-P(Y=1|\bm{x})} = \bm{x}^T\bm{\beta}$ の係数推定量 $\hat{\bm{\beta}}$ は、最尤推定量として漸近正規性を持ちます。
$$ \hat{\bm{\beta}} \overset{\text{approx}}{\sim} N(\bm{\beta}, \bm{I}(\bm{\beta})^{-1}) $$
ここで $\bm{I}(\bm{\beta})$ はフィッシャー情報行列です。
オッズ比 $e^{\beta_j}$ の信頼区間を求めたいとき、$g(\beta_j) = e^{\beta_j}$ に対してデルタ法を適用します。
$$ g'(\beta_j) = e^{\beta_j} $$
$$ \text{Var}(e^{\hat{\beta}_j}) \approx (e^{\beta_j})^2 \cdot \text{Var}(\hat{\beta}_j) $$
ただし実用上は、$\hat{\beta}_j$ に対して信頼区間を構成し、端点に指数関数を適用する方が一般的です。
$$ e^{\hat{\beta}_j \pm z_{\alpha/2}\sqrt{\text{Var}(\hat{\beta}_j)}} $$
Pythonでの実装
実装1：1変数デルタ法の検証
指数分布からの標本平均に対数変換を適用し、デルタ法の近似精度を検証します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# 指数分布 Exp(lambda) の標本平均の対数
# E[X] = 1/lambda, Var(X) = 1/lambda^2
# g(x) = log(x), g'(x) = 1/x
# デルタ法: Var(log(X_bar)) ≈ Var(X_bar) / mu^2 = (1/lambda^2) / (n * (1/lambda)^2) = 1/n
lam = 2.0
mu = 1.0 / lam # 真の期待値
sigma2 = 1.0 / lam**2 # 真の分散
ns = [5, 10, 20, 50, 100, 500]
n_sim = 50000
np.random.seed(42)
fig, axes = plt.subplots(2, 3, figsize=(16, 10))
axes = axes.ravel()
for idx, n in enumerate(ns):
ax = axes[idx]
# シミュレーション: log(X_bar)のサンプリング分布
log_means = np.array([
np.log(np.mean(np.random.exponential(1.0/lam, n)))
for _ in range(n_sim)
])
# デルタ法の予測
# E[log(X_bar)] ≈ log(mu) + bias
mean_delta = np.log(mu)
# Var(log(X_bar)) ≈ sigma^2 / (n * mu^2) = (1/lambda^2) / (n * (1/lambda)^2) = 1/n
var_delta = sigma2 / (n * mu**2)
# 2次デルタ法のバイアス: g''(mu)*sigma^2/(2n) = (-1/mu^2)*(sigma^2)/(2n)
bias_2nd = -sigma2 / (2 * n * mu**2)
# ヒストグラム
ax.hist(log_means, bins=60, density=True, alpha=0.6,
color='steelblue', edgecolor='white', label='Simulation')
# デルタ法による正規近似
x = np.linspace(mean_delta - 4*np.sqrt(var_delta),
mean_delta + 4*np.sqrt(var_delta), 300)
ax.plot(x, stats.norm.pdf(x, mean_delta, np.sqrt(var_delta)),
'r-', linewidth=2.5, label='Delta method (1st)')
# 2次デルタ法（バイアス補正付き）
ax.plot(x, stats.norm.pdf(x, mean_delta + bias_2nd, np.sqrt(var_delta)),
'g--', linewidth=2, label='Delta method (2nd)')
ax.set_title(f'n = {n}', fontsize=13)
ax.set_xlabel(r'$\log(\bar{X})$', fontsize=11)
ax.set_ylabel('Density', fontsize=11)
ax.legend(fontsize=8)
# 統計量の表示
ax.text(0.05, 0.95,
f'Sim mean={np.mean(log_means):.4f}\n'
f'Delta mean={mean_delta:.4f}\n'
f'2nd mean={mean_delta+bias_2nd:.4f}\n'
f'Sim var={np.var(log_means):.5f}\n'
f'Delta var={var_delta:.5f}',
transform=ax.transAxes, fontsize=8, verticalalignment='top',
bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
plt.suptitle(r'Delta Method for $\log(\bar{X})$, $X \sim Exp(2)$',
fontsize=15, y=1.01)
plt.tight_layout()
plt.savefig("delta_method_log.png", dpi=150, bbox_inches="tight")
plt.show()
$n$ が小さいとき（$n=5$ など）は対数変換による歪みが残り、1次デルタ法の正規近似はやや不正確ですが、2次デルタ法の期待値補正により改善されます。$n$ が大きくなるにつれ、デルタ法の近似はほぼ完全に一致します。
実装2：多変数デルタ法（比率の分散）
2つの独立な正規母集団からの標本平均の比 $R = \bar{X}/\bar{Y}$ に対して多変数デルタ法を適用し、シミュレーションと比較します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# パラメータ
mu1, sigma1 = 10.0, 3.0 # X の母集団
mu2, sigma2 = 5.0, 2.0 # Y の母集団
R_true = mu1 / mu2 # 真の比
ns = [10, 30, 100, 500]
n_sim = 50000
np.random.seed(123)
fig, axes = plt.subplots(2, 2, figsize=(12, 9))
axes = axes.ravel()
for idx, n in enumerate(ns):
ax = axes[idx]
# シミュレーション
ratios = np.array([
np.mean(np.random.normal(mu1, sigma1, n)) /
np.mean(np.random.normal(mu2, sigma2, n))
for _ in range(n_sim)
])
# 多変数デルタ法
# Var(R) ≈ (1/mu2^2)(sigma1^2/n) + (mu1^2/mu2^4)(sigma2^2/n)
var_delta = (sigma1**2 / mu2**2 + mu1**2 * sigma2**2 / mu2**4) / n
# ヒストグラム
ax.hist(ratios, bins=80, density=True, alpha=0.6,
color='steelblue', edgecolor='white', label='Simulation')
# デルタ法の正規近似
x = np.linspace(R_true - 4*np.sqrt(var_delta),
R_true + 4*np.sqrt(var_delta), 300)
ax.plot(x, stats.norm.pdf(x, R_true, np.sqrt(var_delta)),
'r-', linewidth=2.5, label='Delta method')
ax.set_title(f'n = {n}', fontsize=13)
ax.set_xlabel(r'$R = \bar{X}/\bar{Y}$', fontsize=11)
ax.set_ylabel('Density', fontsize=11)
ax.legend(fontsize=10)
ax.axvline(R_true, color='green', linestyle='--', alpha=0.7)
# 分散の比較
ax.text(0.05, 0.95,
f'Sim var={np.var(ratios):.5f}\nDelta var={var_delta:.5f}',
transform=ax.transAxes, fontsize=9, verticalalignment='top',
bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
plt.suptitle(r'Multivariate Delta Method: $R = \bar{X}/\bar{Y}$',
fontsize=15, y=1.01)
plt.tight_layout()
plt.savefig("delta_method_ratio.png", dpi=150, bbox_inches="tight")
plt.show()
実装3：オッズ比の信頼区間
$2 \times 2$ 分割表からオッズ比の信頼区間をデルタ法で計算します。
import numpy as np
def odds_ratio_ci(a, b, c, d, alpha=0.05):
"""
2x2分割表からオッズ比の点推定と信頼区間を
デルタ法（対数オッズ比の正規近似）で計算
| Disease+ | Disease- |
Treatment | a | b |
Control | c | d |
"""
# オッズ比の点推定
OR = (a * d) / (b * c)
# 対数オッズ比
log_OR = np.log(OR)
# デルタ法による対数オッズ比の標準誤差
se_log_OR = np.sqrt(1/a + 1/b + 1/c + 1/d)
# 信頼区間（対数スケール）
z = -1 * __import__('scipy').stats.norm.ppf(alpha / 2)
ci_log_lower = log_OR - z * se_log_OR
ci_log_upper = log_OR + z * se_log_OR
# 元のスケールに戻す
ci_lower = np.exp(ci_log_lower)
ci_upper = np.exp(ci_log_upper)
return OR, (ci_lower, ci_upper), se_log_OR
# 例: 臨床試験データ
# 処置群: 30人中20人が回復
# 対照群: 30人中10人が回復
a, b = 20, 10 # 処置群: 回復, 非回復
c, d = 10, 20 # 対照群: 回復, 非回復
OR, ci, se = odds_ratio_ci(a, b, c, d)
print("=== オッズ比の推定（デルタ法） ===")
print(f"2x2 分割表:")
print(f" 処置群: 回復={a}, 非回復={b}")
print(f" 対照群: 回復={c}, 非回復={d}")
print(f"\nオッズ比 OR = {OR:.4f}")
print(f"log(OR) = {np.log(OR):.4f}")
print(f"SE[log(OR)] = {se:.4f}")
print(f"95% 信頼区間: ({ci[0]:.4f}, {ci[1]:.4f})")
print(f"\nOR > 1 かつ信頼区間が1を含まない → 処置効果あり" if ci[0] > 1 else
f"\n信頼区間が1を含む → 有意な処置効果なし")
実装4：デルタ法 vs ブートストラップ法の比較
デルタ法による分散推定の精度を、ブートストラップ法と比較検証します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# 設定: ポアソン分布からの標本平均の平方根の分散
# X ~ Poi(lambda), E[X] = lambda, Var(X) = lambda
# g(x) = sqrt(x), g'(x) = 1/(2*sqrt(x))
# デルタ法: Var(g(X_bar)) ≈ [g'(lambda)]^2 * lambda/n = 1/(4n)
lam = 4.0
ns = [10, 20, 50, 100, 200, 500]
n_sim = 5000
n_boot = 2000
np.random.seed(456)
# 各方法による標準誤差の推定を比較
results = {
'n': [],
'true_se': [], # シミュレーションで求めた真のSE
'delta_se': [], # デルタ法のSE
'boot_se_mean': [], # ブートストラップSEの平均
'boot_se_std': [], # ブートストラップSEの標準偏差
}
for n in ns:
# 真のSE（大量シミュレーション）
g_values = np.array([
np.sqrt(np.mean(np.random.poisson(lam, n)))
for _ in range(50000)
])
true_se = np.std(g_values)
# デルタ法のSE
# Var(sqrt(X_bar)) ≈ [1/(2*sqrt(lambda))]^2 * lambda/n = 1/(4n)
delta_se = np.sqrt(1.0 / (4 * n))
# ブートストラップSEの分布
boot_ses = []
for _ in range(n_sim):
sample = np.random.poisson(lam, n)
# ブートストラップ
boot_stats = np.array([
np.sqrt(np.mean(np.random.choice(sample, n, replace=True)))
for _ in range(n_boot)
])
boot_ses.append(np.std(boot_stats))
results['n'].append(n)
results['true_se'].append(true_se)
results['delta_se'].append(delta_se)
results['boot_se_mean'].append(np.mean(boot_ses))
results['boot_se_std'].append(np.std(boot_ses))
# 結果の表示
print("=" * 70)
print(f"{'n':>6} | {'True SE':>10} | {'Delta SE':>10} | {'Boot SE (mean)':>15} | {'Boot SE (std)':>14}")
print("-" * 70)
for i in range(len(ns)):
print(f"{results['n'][i]:>6} | {results['true_se'][i]:>10.5f} | "
f"{results['delta_se'][i]:>10.5f} | {results['boot_se_mean'][i]:>15.5f} | "
f"{results['boot_se_std'][i]:>14.5f}")
# 可視化
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))
# SEの比較
ax1.plot(ns, results['true_se'], 'ko-', linewidth=2, markersize=8, label='True SE')
ax1.plot(ns, results['delta_se'], 'r^--', linewidth=2, markersize=8, label='Delta method')
ax1.errorbar(ns, results['boot_se_mean'], yerr=results['boot_se_std'],
fmt='bs-.', linewidth=2, markersize=8, capsize=5, label='Bootstrap')
ax1.set_xlabel('Sample size n', fontsize=12)
ax1.set_ylabel('Standard Error', fontsize=12)
ax1.set_title(r'SE of $\sqrt{\bar{X}}$ where $X \sim Poi(4)$', fontsize=14)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xscale('log')
# 相対誤差
delta_rel_err = np.abs(np.array(results['delta_se']) - np.array(results['true_se'])) / np.array(results['true_se'])
boot_rel_err = np.abs(np.array(results['boot_se_mean']) - np.array(results['true_se'])) / np.array(results['true_se'])
ax2.plot(ns, delta_rel_err * 100, 'r^--', linewidth=2, markersize=8, label='Delta method')
ax2.plot(ns, boot_rel_err * 100, 'bs-.', linewidth=2, markersize=8, label='Bootstrap')
ax2.set_xlabel('Sample size n', fontsize=12)
ax2.set_ylabel('Relative Error (%)', fontsize=12)
ax2.set_title('Relative Error of SE Estimation', fontsize=14)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_xscale('log')
plt.tight_layout()
plt.savefig("delta_vs_bootstrap.png", dpi=150, bbox_inches="tight")
plt.show()
結果から以下のことが読み取れます。
- デルタ法 は閉じた形の公式を与え、計算コストが極めて低いです。$n$ が中程度以上（$n \geq 30$ 程度）であれば、十分な精度が得られます。
- ブートストラップ法 はデルタ法と同程度の精度を達成しますが、計算コストが大きく、また推定値自体にばらつき（ブートストラップ標準誤差の標準偏差）があります。
- $n$ が小さいとき、デルタ法は線形近似の誤差が大きくなり、ブートストラップ法の方が信頼性が高い場合があります。
実装5：信頼区間の被覆確率
デルタ法による信頼区間が名目上の被覆確率を達成しているかを検証します。
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# 指数分布の標本平均の逆数 g(X_bar) = 1/X_bar
# E[X] = 1/lambda, Var(X) = 1/lambda^2
# g'(x) = -1/x^2
# Var(g(X_bar)) ≈ (1/mu^2)^2 * sigma^2/n = lambda^4 * (1/lambda^2)/n = lambda^2/n
lam = 3.0
mu = 1.0 / lam
sigma2 = 1.0 / lam**2
ns = [10, 20, 30, 50, 100, 200, 500, 1000]
n_sim = 10000
alpha = 0.05
z = stats.norm.ppf(1 - alpha/2)
np.random.seed(789)
coverage_delta = []
for n in ns:
cover_count = 0
for _ in range(n_sim):
sample = np.random.exponential(1.0/lam, n)
x_bar = np.mean(sample)
g_hat = 1.0 / x_bar # g(X_bar) = 1/X_bar
# デルタ法のSE: |g'(X_bar)| * SE(X_bar)
# g'(x) = -1/x^2, SE(X_bar) = s/sqrt(n)
s = np.std(sample, ddof=1)
se_delta = (1.0 / x_bar**2) * (s / np.sqrt(n))
# 信頼区間
ci_lower = g_hat - z * se_delta
ci_upper = g_hat + z * se_delta
# 真の値 lambda = 3.0 が信頼区間に入っているか
if ci_lower <= lam <= ci_upper:
cover_count += 1
coverage_delta.append(cover_count / n_sim)
# 可視化
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(ns, coverage_delta, 'bo-', linewidth=2, markersize=8, label='Delta method CI')
ax.axhline(y=1-alpha, color='red', linestyle='--', linewidth=2,
label=f'Nominal level ({1-alpha:.2f})')
ax.fill_between(ns, 1-alpha-0.02, 1-alpha+0.02, alpha=0.1, color='red')
ax.set_xlabel('Sample size n', fontsize=12)
ax.set_ylabel('Coverage Probability', fontsize=12)
ax.set_title(r'Coverage of Delta Method CI for $1/\bar{X}$, $X \sim Exp(3)$',
fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
ax.set_ylim(0.85, 1.0)
ax.set_xscale('log')
plt.tight_layout()
plt.savefig("delta_coverage.png", dpi=150, bbox_inches="tight")
plt.show()
print("被覆確率:")
for n, cov in zip(ns, coverage_delta):
print(f" n={n:5d}: {cov:.4f}")
$n$ が10程度の小標本ではデルタ法の信頼区間の被覆確率は名目水準 95% を下回りますが、$n$ が50を超えるあたりから95%に近づき、$n \geq 100$ ではほぼ正確に名目水準を達成することが確認できます。
まとめ
本記事では、デルタ法の理論と応用について包括的に解説しました。
- デルタ法の動機として、漸近正規な推定量の関数の分散を求めるという問題を設定しました
- 1変数デルタ法をテイラー展開とスラツキーの定理で厳密に導出し、$\text{Var}(g(T_n)) \approx [g'(\theta)]^2\sigma^2/n$ を得ました
- 多変数デルタ法を導出し、$\text{Var}(g(\bm{T}_n)) \approx \nabla g^T \bm{\Sigma}\,\nabla g / n$ を示しました
- 2次デルタ法によりバイアス $\frac{g”(\theta)\sigma^2}{2n}$ の近似式を導出しました
- 対数変換、比率の分散、オッズ比の信頼区間、ロジスティック回帰への具体的応用を示しました
- Pythonでデルタ法の近似精度を検証し、ブートストラップ法との比較で計算コストと精度のトレードオフを明らかにしました
- 信頼区間の被覆確率シミュレーションにより、$n \geq 50$ 程度でデルタ法の近似が実用的に十分であることを確認しました
デルタ法は統計学のあらゆる場面で登場する基本ツールです。次のステップとして、以下の記事も参考にしてください。