はじめに
東京大学/株式会社Nospare リサーチャーの栗栖です．
この記事から数回に分けて，経験尤度法とその多様体データ解析への応用について紹介します．今回の記事では経験尤度法の考え方について紹介します．
経験尤度法
経験尤度法 (empirical likelihood method, EL法)は統計学者 Art, B. Owen によって提案された方法(Owen (1988, 1990, 1991))で，現在では統計学を含め，計量経済学などの分野でも用いられています．EL法はモーメント制約モデルにおけるパラメータの推定，推測のための方法で，その名前の通り最尤推定法と密接に関係しています．最尤法はデータの同時確率密度を最大化するようにパラメータの値を決める推定手法であり，その実行には特定のパラメータにより決定されるパラメトリックな分布を仮定する必要がありました．経験尤度法は分布の形に特定の仮定を置かないノンパラメトリックな最尤法(nonparametric MLE, NPMLE)と理解することができます．
ノンパラメトリック最尤推定量
まずノンパラメトリックな尤度の定義から始めます．$W_1,\dots, W_n$をi.i.d.確率ベクトルとし，それぞれの値が実現する確率を$p_1,\dots,p_n$であるような分布を考えます．このとき，$W_1,\dots,W_n$の尤度関数は同時確率$\prod_{i=1}^{n}p_i$であり，対数尤度は$\sum_{i=1}^{n}\log p_i$となります．この値を最大化するような$p_i$がノンパラメトリックMLEです．ただし，確率は非負でなければならず，それらの和は1でなければならないので以下のような制約付き最大化問題を考えます：
\begin{align*}
\max_{p_1,\dots,p_n}\sum_{i=1}^{n}\log p_i\quad \text{s.t.}\ \sum_{i=1}^{n}p_i=1,\ p_i \geq 0,\ i=1,\dots,n.
\end{align*}
Lagrange の未定乗数法を用いて上記の最適化問題を解いてみます．以下の関数を考えましょう：
f(p_1,\dots,p_n, \kappa) = \sum_{i=1}^{n} \log p_i - \kappa \left(\sum_{i=1}^{n}p_i - 1\right)
ここで，$\lim_{x \to 0}\log x = -\infty$となるので，上記の右辺第1項は$p_i \geq 0$の制約も含むことに注意しておきます．この最適化問題の解は以下の1階条件を満たします：
\begin{align*}
{\partial \over \partial \kappa}f(p_1,\dots,p_n, \kappa) &= -\left(\sum_{i=1}^{n}p_i - 1\right) = 0\\
{\partial \over \partial p_j}f(p_1,\dots,p_n, \kappa) &= {1 \over p_j} - \kappa =0,\ j=1,\dots,n
\end{align*}
これらを解くと，$\kappa = n$, $\hat{p}_i = 1/n$, $i=1,\dots,n$となり，各実現値に$1/n$の確率を割り当てたものが確率関数のノンパラメトリックMLEということになります．
EL推定量
次にノンパラメトリック尤度の構成にデータが満たすべきモーメント条件を組み込むことを考えましょう．母集団の平均や分散，線形回帰モデルの回帰係数など，興味のある真のパラメータ$\theta_0 \in \Theta
\subset \mathbb{R}^p$とし，そのパラメータがモーメント条件$E_{\theta_0}[g(W_1;\theta_0)] = 0 \in \mathbb{R}^q$を満たすとします．ここで，$E_{\theta_0}$は真のパラメータの下での期待値を意味します．
モーメント条件の例
(1) 平均の推定に興味がある場合：$g(w;\theta) = w-\theta$として
E_{\theta_0}[g(W_1;\theta_0)]= E_{\theta_0}[W_1] - \theta_0 = 0.
(2) 平均($\mu$)と分散($\sigma^2$)の推定に興味がある場合：$\theta_0 = (\mu,\sigma^2)'$, $g(w;\theta) = (w-\theta, (w - \mu)^2 - \sigma^2)'$として
E_{\theta_0}[g(W_1;\theta_0)]= (E_{\theta_0}[W_1] - \mu, E_{\theta_0}[(W_1-\mu)^2]-\sigma^2)' = (0,0)'.
(3) 線形回帰モデルの回帰係数の推定に興味がある場合：$W_i = (Y_i,X_i)' \in \mathbb{R} \times \mathbb{R}^p$, $Y_i = X_i'\beta + \varepsilon_i$, $E[X_i\varepsilon_i]=0$とすると，$\theta_0 = \beta \in \mathbb{R}^p$として
E_{\theta_0}[g(W_1;\theta_0)]= E[X_i(Y_i - X'_i\beta)]=0 \in \mathbb{R}^p.
EL推定量の構成
モーメント条件をノンパラメトリックな尤度最大化問題に組み込むこと考えましょう．データ$W_i$が実現する確率が$p_i$の場合，モーメント条件$E_{\theta_0}[g(W_1;\theta_0)] = 0$の標本対応は
\sum_{i=1}^{n}p_i g(W_i;\theta_0) = 0
と表されます．$\theta$を固定した下で，上記のモーメント条件を制約とした尤度最大化問題を解き，その最大値を$\ell(\theta)$とします．即ち，
\begin{align*}
\ell(\theta):=\max_{p_1,\dots,p_n}\sum_{i=1}^{n}\log p_i \quad \text{s.t.}\ \sum_{i=1}^{n}p_i=1,\ \sum_{i=1}^{n}p_i g(W_i;\theta) = 0.
\end{align*}
$\ell(\theta)$は$\theta$の関数であり，この関数を$\theta$に関して最大化する$\hat{\theta}_{EL}$が$\theta_0$のEL推定量となります．Lagrange の未定乗数法を用いて上記の最適化問題を解いてみます．以下の関数を考えましょう：
f(p_1,\dots,p_n, \kappa) = \sum_{i=1}^{n} \log p_i - \kappa \left(\sum_{i=1}^{n}p_i - 1\right) - n\lambda'\sum_{i=1}^{n}p_i g(W_i;\theta)
この最適化問題の1階条件は以下で与えられます：
\begin{align}
{\partial \over \partial p_j}f(p_1,\dots,p_n, \kappa, \lambda) &= {1 \over p_j} - \kappa -n\lambda'g(W_i;\theta) =0,\ j=1,\dots,n \label{FOC-EL1}\tag{1} \\
{\partial \over \partial \kappa}f(p_1,\dots,p_n, \kappa, \lambda) &= -\left(\sum_{i=1}^{n}p_i - 1\right) = 0 \label{FOC-EL2}\tag{2}\\
{\partial \over \partial \lambda}f(p_1,\dots,p_n, \kappa, \lambda) &= -n\sum_{i=1}^{n}p_ig(W_i;\theta) = 0 \label{FOC-EL3}\tag{3}
\end{align}
以上の条件より，
\begin{align*}
(\ref{FOC-EL1}) &\Leftrightarrow 1-\kappa p_i - n\lambda' p_i g(W_i;\theta) = 0 \Rightarrow \sum_{i=1}^{n}\left(1-\kappa p_i - n\lambda' p_i g(W_i;\theta)\right) = 0\\
&\stackrel{(\ref{FOC-EL2})}{\Rightarrow} n - \kappa - n\lambda' \sum_{i=1}^{n}p_i g(W_i;\theta) = 0 \stackrel{(\ref{FOC-EL3})}{\Rightarrow} n - \kappa = 0.
\end{align*}
となるので，
p_i = {1 \over n}{1 \over 1 + \lambda' g(W_i;\theta)},\ i=1,\dots,n.
これらを($\ref{FOC-EL3}$)に代入すると，$\lambda$は
g(\lambda) = {1 \over n}\sum_{i=1}^{n}{g(W_i;\theta) \over 1 + \lambda' g(W_i;\theta)} = 0
の解であり，その解を$\lambda(\theta)$として改めて$p_i$に代入すると，
p_i = {1 \over n}{1 \over 1 + \lambda(\theta)' g(W_i;\theta)},\ i=1,\dots,n.
以上の議論により，
\ell(\theta) = -n\log n - \sum_{i=1}^{n}\log(1 + \lambda(\theta)' g(W_i;\theta))
と表されます．第1項は$\theta$に依存しないことから，$\ell(\theta)$から第1項を除いた関数を$Q(\theta)$とすると，EL推定量は
\begin{align*}
\hat{\theta}_{EL} &= \mathrm{argmax}_{\theta \in \Theta}\left\{-\sum_{i=1}^{n}\log (1 + \lambda(\theta)'g(W_i;\theta))\right\}\\
&= \mathrm{argmax}_{\theta \in \Theta}Q(\theta)
\end{align*}
の形で与えられます．またEL推定量を$p_i$に代入すると
\hat{p}_i = {1 \over n}{1 \over 1 + \lambda(\hat{\theta}_{EL})' g(W_i; \hat{\theta}_{EL})},\ i=1,\dots,n.
が求められ，$W_1,\dots, W_n$の分布を推定することもできます．
EL推定量の性質
モーメント条件$E_{\theta_0}[g(W_1;\theta_0)] = 0$が正しいと仮定します．このとき，EL推定量の定義に現れる関数$Q$に対して，適当な条件の下で以下の結果が成り立つことが知られています：
\begin{align*}
Q(\theta_0) &\stackrel{d}{\to} \chi_p^2.
\end{align*}
ここで，$\chi_p^2$は自由度$p$のカイ2乗分布です．従って，$q_{1-\alpha}$を$\chi_p^2$分布の$1-\alpha$分位点とすると，上記の結果を用いて$Q(\theta)$のレベル集合 $\{\theta \in \Theta: Q(\theta) \leq q_{1-\alpha}\}$により$\theta_0$の$100(1-\alpha)$%信頼領域を構成することができます．また$q_{1-\alpha}$を critical value とすることで，帰無仮説$\mathbb{H}_0: \theta_0 = \theta$ (対立仮説$\mathbb{H}_1: \theta_0 \neq \theta$) の有意水準$100\alpha$%の検定を実行することができます．
経験尤度法のメリット・デメリット
最後に経験尤度法のメリット・デメリットとして主要なものをまとめておきます．
メリット：
経験尤度法では，推定量の漸近正規性に基づく信頼領域の構成や仮説検定と異なり，信頼領域の構成や仮説検定において$Q(\theta_0)$の漸近分散(または局外母数)の推定が必要ない．この記事では詳細は省きますが，EL推定量は適当な条件の下で，一般化モーメント法(generalized method of moments, GMM)などの方法と比較して理論的に良い性質をもつことが知られています．
デメリット：
経験尤度法ではEL推定量を最小化問題の解として数値的に求める必要がある．
まとめ
この記事では経験尤度法について紹介しました．次回の記事では経験尤度法を用いて多様体データの分析を行う方法(Kurisu and Otsu(2024))について紹介します．株式会社Nospareには統計学の様々な分野を専門とする研究者が所属しています．統計アドバイザリーやビジネスデータ分析につきましては株式会社Nospareまでお問い合わせください．
参考文献
[1] Kurisu, D. and Otsu, T. (2024) Empirical likelihood for manifolds.
[2] Owen, A. B. (1988) Empirical likelihood ratio confidence intervals for a single functional. Biometrika 75, 237-249.
[3]Owen, A. B. (1990) Empirical likelihood confidence regions. Annals of Statistics 18, 90-120.
[4]Owen, A. B. (1991) Empirical likelihood for linear models. Annals of Statistics 19, 1725-1747.