こんにちは，株式会社Nospare・明治大学の小林です．本記事では尤度関数が得られない場合に事後分布を近似的に求める方法について簡単に紹介します．
尤度関数フリーの手法
ベイズ分析で関心の対象であるパラメータ$\theta$の事後分布は
p(\theta|y)=\frac{f(y|\theta)p(\theta)}{f(y)}\propto f(y|\theta)p(\theta)
で表記されます（ここで尤度関数$f(y|\theta)$，事前分布$p(\theta)$，データの周辺分布$f(y)$）．標準的な多くの統計モデルでは尤度関数$f(y|\theta)$は解析的に得られたり，その数値的な評価は簡単に行うことができます．例えば，$p(\theta|y)$に対してメロポリス・ヘイスティングス（MH）アルゴリズムを適用しようとすると，採択確率を計算する際には次のMH比
\frac{p(\theta^*|y)q(\theta|\theta^*)}{p(\theta|y)q(\theta^*|\theta)}=\frac{f(y|\theta^*)p(\theta^*)q(\theta|\theta^*)}{f(y|\theta)p(\theta)q(\theta^*|\theta)}
を評価します（$q(\theta^* | \theta)$は提案分布，$\theta^*$は提案値）が，これには尤度関数，事前分布$p(\theta)$，提案分布の評価が必要になります．
一方で，例えば微分方程式に基づくモデル，一部の空間モデル，分位点や特性関数によって定義される確率分布などでは，モデルに対応する尤度関数が解析的に得られなかったり尤度関数の評価コストが非常に高くなったりすることがあります．このようなモデルは疫学・化学・環境学など様々な分野において頻繁に用いられます．尤度関数が評価できない場合には，上のMH比を評価することはできませんし，評価することができたとしても計算コストが非常に高い場合には現実的な時間内でモデルの推定が不可能になってしまいます．
簡単な例として，$g$-and-$k$分布（$g$-and-$h$もあります）を紹介します．$g$-and-$k$分布は次の分位点によって定義されます：
Q(u)=A+B(1+c \tanh(g z_u/2))z_u(1+z_u^2)^k,\quad u\in(0,1)
ここで$z_u\Phi(u)^{-1}$は標準正規分布の分位点，$A$, $B$, $c$, $g$, $k$はパラメータです．$f(y)=\frac{1}{q(Q^{-1}(y))}$のように分位点から確率密度関数を導出することも可能ですが（$q(u)=dQ(u)/du$），解析的な形は得られません（数値的には評価することは可能です）．
尤度関数フリー（likelihood-free; LF）の手法は，このような状況において近似的にベイズ推定を行うための手法の広いクラスであり，ここ20年近くで多くの研究の蓄積がされてきています．Pacchiardi and Dutta (2022)は，尤度関数フリーの方法の紹介においておおまかに2つのジャンルに分けています．
- 近似ベイズ計算（ABC; approximate Bayesian computation）
- サロゲート尤度
サロゲート尤度の方には合成尤度（SL; synthetic likelihood）と密度比（DR; density ratio）推定の方法が含まれます．これらの方法を適用するには，確率モデルからのデータを生成することができる，ということが必要になります．上述に挙げたようなモデルたちや，他のコンピュータモデル等は尤度関数が得られませんがモデルからのシミュレーションは可能なので，尤度関数フリーの手法が適用できます．上で紹介した$g$-and-$k$分布からのシミュレーションは簡単に行うことができ，一様分布からのシミュレーション$u\sim U(0,1)$を分位点に代入するだけです．
近似ベイズ計算
ABCは尤度フリーの方法で一番標準的な方法といえます．ABCは，観測されたデータ$y$と近いようなデータを生成することができるようなパラメータを事後分布からのサンプルとしてみなします．とあるパラメータの値のもとで人工データ$x$を生成し，$y$と$x$の近さは，通常それぞれの要約統計量$s(y)$と$s(x)$の間の距離$d(s(y),s(x))$で測り，その距離が一定の値$\epsilon$を下回る場合にそのパラメータの値は事後分布からのサンプルとして保持します．このABCのもとでの近似事後分布は
p_\epsilon(\theta|s(y))\propto \pi(\theta)\int I\left\{d(s(y),s(x))<\epsilon\right\}f(x|\theta)dx
で与えられます．
ABCを適用するには
- 要約統計量$s$
- 距離関数$d$
- 閾値$\epsilon$
を選択する必要があり，選択の方法に関して様々なものが提案されてきています．またアルゴリズムに関しても最も単純な棄却サンプリングからMHアルゴリズム，逐次モンテカルロ法などに基づくものが提案されてきています．
サロゲート尤度
サロゲート尤度のアプローチは何らかの形で近似尤度を構成します．
合成尤度
合成尤度は評価できない尤度関数を要約統計量に対する正規分布で置き換えます：
s(y)|\theta \sim N(\mu_\theta,\Sigma_\theta)
ここで$\mu_\theta$と$\Sigma_\theta$は，パラメータの値を所与にモデルから生成されたデータ$x$に基づいて計算されます．この近似的な正規尤度を上のMH比の$f(y| \theta)$へ代入することで近似的な事後分布からのサンプリングを行うことができます．
ABCでは一つのパラメータの値に対して$x$の生成は1回で十分ですが，合成尤度の場合は正規分布の平均と分散を計算するために$x$の生成を複数回行う必要があります．要約統計量が正規分布に従うという仮定が制約的であると考えられる場合には，$t$分布や非対称な分布に置き換えられることがあります．
密度比推定
密度比推定のアプローチでは尤度関数とデータの周辺分布の比を推定します：
r(y,\theta)=\frac{f(y|\theta)}{f(y)}
$r(y,\theta)\propto f(y| \theta)$であるので，
\frac{r(y,\theta^*)q(\theta|\theta^*)}{r(y,\theta)q(\theta^*|\theta)}
をMH比として用いることができます．密度比$r(y,\theta)$を推定するためには，$f(y|\theta)$と$f(y)$からのデータの生成を行います．そしてロジスティック回帰によって2種類の人工データを分類し，$r(y,\theta)\approx \exp(h(y,\theta))$とします（$h$はロジスティック回帰の回帰関数）．この方法による尤度関数の構成は合成尤度よりも一般的であり，最近では関連する手法もさらにいくつか提案されてきています．
未知パラメータに依存する正規化定数
ここまでは尤度関数全体が得られない一般的な場合について取り上げてきましたが，関連する特別な状況として未知の正規化定数が未知パラメータに依存する場合（doubly intractable）が挙げられます．この場合でも基本的通常のMCMCなどが適用できないため注意が必要です．詳細は評価できない正規化定数がある場合のベイズ推定について【１】で紹介していますが，上述の尤度関数フリーの方法もこの場合に適用することは可能です．
まとめ
分野によっては尤度関数が得られない確率モデルを扱う状況にはなかなか直面しないかもしれませんが，知っておくと便利だと思われます．本記事で紹介した手法で鍵となるのはモデルからのシミュレーションが可能であることです．doubly intracable固有の方法を適用する際にはシミュレーションは必ずしも必要ではありません．また要約統計量の選択はABCやSLのパフォーマンスに大きな影響を与えるので，ドメイン知識に基づいて注意深く選ぶか，ABCの文脈で提案されている半自動的に要約統計量を構成する手法を使用することが勧められます．
一緒にお仕事しませんか！
株式会社Nospareではベイズ統計学に限らず統計学の様々な分野を専門とする研究者が所属しており，新たな知見を日々追求しています．統計アドバイザリーやビジネスデータの分析につきましては弊社までお問い合わせください．インターンや正社員も随時募集しています！
