パラメータに制約があるベイズ推論

# パラメータ制約
こんにちは，[株式会社Nospare](https://bit.ly/38WtGZH)・明治大学の小林です．
本記事ではモデルパラメータに制約が課された状況でのベイズ推論の方法に関する最近の研究を2つ紹介します．パラメータに対する制約が必要な状況というのは以外にも多く存在します．例えば

* 回帰関数に対する単調制約
* 行列の特定の性質を満たすような制約（正定値性，直交性など）
* 経済学などいった背後にある理論に基づく制約

などが挙げられ，パラメータ$\theta\in R$に対する制約の例としては

* $\theta>0$
* $||\theta||^2=1$（球上）
* $\theta^\top\theta=I$（直交）
* $\theta^\top b = a$

などといったような制約が課されます．制約が課された場合のパラメータ空間$D$は制約なしのパラメータ空間のに含まれる形となります（$D\subset R)$）

通常の制約なしの事後分布が

```math
\pi_R (\theta|Y)\propto \pi_R(\theta)L(\theta;Y)
```

で表されるのに対し，制約付きの事後分布は

```math
\pi_D (\theta|Y)\propto \pi_D(\theta)L(\theta;Y)
```

で表されます．

事後分布は$\theta$に関する不確実性の定量化に必要なものですが，制約を課することによって以下のような問題が生じます：

* 制約を満たすような事前分布のモデリングが難しく，特に$D$がユークリッド空間ではない場合には$D$上の確率分布の選択肢がかなり限定される（例えば$D$が球上であれば殆どの場合はvon Mises分布に関連する分布）
* MCMCなどの数値計算がうまく機能しない

このような問題を解決するためのアプローチはいままでも研究がなされてきてはいますが，本記事では比較的最近提案された方法を2つ

* 制約を緩める（constraint relaxation）
* 事後分布を射影する（posterior projection）

紹介します．

# 制約を緩める（CR; constraint relaxation）
この方法はDuan etal. (2020)[[paper](https://academic.oup.com/biomet/article/107/1/191/5686747),[arXiv](https://arxiv.org/abs/1801.01525)]によって提案されたものです．上述の制約付き事後分布は，$D$についての指示関数$I_D(\theta)$を使って

```math
\pi_D (\theta|Y)\propto \pi_R(\theta)I_D(\theta)L(\theta;Y)
```

のように書くことができると仮定します．この場合，事前分布は通常の制約なしの場合の事前分布$\pi_R(\theta)$に指示関数がかかったものになります．CRでは，この指示関数によるシャープな制約を次のように緩い制約に置き換えます．

```math
\tilde{\pi}_\lambda (\theta|Y)\propto L(\theta;Y)\pi_R(\theta)\exp\left\{-\lambda^{-1}||\nu_D(\theta)||\right\}
```

ここで$||\nu_D(\theta)||$は$\theta$と制約された空間との距離を表します．事前分布$\pi_R(\theta)$は正則で，$D$との空でない共通部分をもつ$R$上の集合上にサポートをもつものとします．

ハイパーパラメータ$\lambda>0$は，事前分布が$D$の周辺にどれだけ集中しているかをコントロールします．$\lambda\rightarrow0$のとき，$\theta\notin D$について$\tilde{\pi}_\lambda (\theta|Y) \rightarrow 0$となり，シャープな制約に近づいていきます．

一方で$\lambda>0$の場合には緩い事後分布$\tilde{\pi}_\lambda (\theta|Y)$は$D$の外にサポートを持ちます．


## 緩い制約の入れ方
$\nu_D(\theta)$のデザイン方法はどのような制約を課したいかに依存します．

* シンプレックスへの制約：$R=[0,1]^r$，$dim(R)=r$，$dim(D)=r-1$，$\nu(\theta)=\sum\theta-1$
* 直線への制約：$R=\mathbb{R}^r$，$dim(R)=r$，$dim(D)=1$，$\nu(\theta)=\theta^\top b$
* 球上への制約：$R=[-1,1]^r$，$dim(R)=r$，$dim(D)=r-1$，$\nu(\theta)=||\theta||^2-1$
* 直交性（Stiefel多様体上へ）の制約：$R=[-1,1]^{n\times k}$，$dim(R)=nk$，$dim(D)=nk-_{k+1}C_2$，$\nu(\theta)=\theta^\top\theta-I$

緩い事後分布のもとでの関数$g(\theta)$の事後期待値はシャープな事後分布のもとでの事後期待値に$\lambda\rightarrow0^+$のときに収束します．またそのオーダーはある条件のもとで$O(\lambda|\log\lambda|)$となります．


このアプローチの利点は，課したい制約に対応する$\nu_D(\theta)$と$\lambda$の値を決めてしまえば，通常の$R$上でのMCMCを適用でき，複雑な問題に対しても適用できるというところにあります．実際，Duan etal. (2020)では，ハミルトニアン・モンテカルロ法を適用しています．またMAP推定値などもアニーリング法などを使って簡単に得ることができます．



# 事後分布を射影する（PP; posterior projection）
こちらの方法は[Sen etal. (2022)](https://arxiv.org/abs/1812.05741)によって提案されたものになります．制約を課さない事後分布$\pi_R(\theta|Y)$を得たあとに，それを$D$上に射影するというのが基本的な考え方になります．$\theta$と$D$との距離を$dist(\theta,D)=\inf({ ||\theta-\tilde{\theta}||:\tilde{\theta}\in D })$で表すとすると，$||\theta-\tilde{\theta}||=dist(\theta,D)$を満たす$D$上の点は，$\theta$を最もよく近似します．これを用いて次の$R$から$D$への射影を定義すると

```math
T\theta = \left\{\tilde{\theta}\in D:||\theta-\tilde{\theta}||=dist(\theta,D) \right\}
```

$B\in D$制約付き事後分布は

```math
\pi_D(\tilde{\theta}|Y)=\pi_R(T^{-1}\tilde{\theta}|Y)
```

で与えられます．Sen etal.(2022)では射影された事後分布の漸近正規性などの理論的な性質を示しています．

この方法の利点としてやはりシャープな制約が保持されることが挙げられ，CRと比べてパラメータを正確に推定することができます．また，制約なしの事後分布からのサンプルに対して射影を行うので，通常の数値計算の手法を適用することができます．


## 正規分布の例
次の簡単な例をとりあげます．

* $Y=(Y_1,\dots,Y_n)$ 
* $Y_i\sim N(\theta,1)$
* $\theta\geq 0$（非負制約）


 まずは$\theta$に対し，非負制約を満たすように切断正規事前分布を仮定します：$\theta\sim N_{[0,\infty)}(0,10^3)$．この切断正規事前分布のもとでは，事後分布は切断正規分布$N_{[0,\infty)}(\theta_n,\sigma_n^2)$となります（$\theta_n=n\bar{Y}(10^{-3}+1)^{-1}）$, $\sigma_n^2=(10^{-3}+1)^{-1}$)．この事後分布では$\theta=0$の確率測度はゼロで，$\theta$の事後期待値は

```math
E[\theta|Y]=\theta_n+\phi(\alpha)+\frac{\sigma_n}{1-\Phi(\alpha)}
```

となります．ここで$\alpha=-\theta_n/\sigma_n$，$\phi$と$\Phi$はそれぞれ標準正規分布の密度関数と分布関数です．

一方でPPのもとでは，まず制約なしの事前分布$\theta\sim N(0,10^3)$を仮定して事後分布を得ます．そして事後分布の負の領域にある$\theta$を制約の境界である$0$に射影します．これにより，射影後の事後分布は

```math
\pi_D(\tilde{\theta}|Y)=\Phi(\alpha)\delta_{\tilde{\theta}=0} + (1-\Phi(\alpha))N_{[0,\infty)}(\theta_n,\sigma_n^2)
```

で与えられます（$\delta_{\tilde{\theta}=0}$は$\tilde{\theta}=0$でのディラック測度）．事後期待値は

```math
E[\tilde{\theta}|Y]=(1-\Phi(\alpha))\left[\theta_n+\phi(\alpha)\frac{\sigma_n}{1-\Phi(\alpha)}\right]
```

となり，$E[\theta|Y]$よりも$\theta_n$に近くなります．


# おわりに
本記事では，パラメータ制約を（近似的にも）事後分布全体が満たすように事後分布を構成する方法を紹介しました．ところが応用問題によっては，分布全体が制約を満たさなくてもモーメントなど，分布のある側面だけが特定の制約を満たしていれば十分な場合もあります．今後はそのような方法についても取り上げたいと思います．

# 参考文献
1. Duan, L., Young, A.L., Nishimura, A. and Dunson, D.B. (2020). Bayesian constraint relaxation. Biometrika, 1, 191-204.
1. Sen, D., Patra, S. and Dunson, D.B. (2022). Constrained inference through posterior projections. arXiv:1812.05741.

__一緒にお仕事しませんか！__
[株式会社Nospare](https://bit.ly/38WtGZH)ではベイズ統計学に限らず統計学の様々な分野を専門とする研究者が所属しており，新たな知見を日々追求しています．ビジネスデータの分析や研究者による統計学関連のセミナーにつきましては弊社までお問い合わせください．インターンや正社員も随時募集しています！


<a href="https://bit.ly/38WtGZH">
<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1023506/d5c625c0-d6f4-680b-d702-3c40f11d6c39.jpeg" alt="株式会社Nospare" width=400>
</a>