切断正規分布に従うパラメータの再パラメータ化

1/20 2025年2月23日 ベイズ統計学勉強会 春'24（ベイズ塾春合宿） 切断正規分布に従うパラメータの 再パラメータ化 武藤 拓之 (Hiroyuki Muto) 大阪公立大学大学院現代システム科学研究科

2/20 再パラメータ化の 復習

3/20 Stanにおける再パラメータ化の必要性 ◼ 再パラメータ化 (reparameterization) ⚫ モデル式に含まれるパラメータを工夫して変数変換することによって， サンプリングを効率的にする方法 (松浦, 2016) ⚫ Stanでモデルを書くときには， 標準正規分布や[0,1]の範囲の一様分布を用いて 再パラメータ化すると効率が改善しやすい。 (ただしモデルによっては悪化することもある。) ⚫ MCMCが全然収束しないときや， divergent transitionsが発生するときは， 再パラメータ化すると解決する場合がある。 (経験上，特に時系列モデルでは効率が劇的に変化することが多い。) 再パラメータ化を自然にできるようになるとつまづきを回避しやすい

4/20 例1）ニールの漏斗 (Neal’s Funnel) ※松浦 (2016, Ch. 9) を参考にした。 ◼ うまくいかないモデル ⚫ 𝑌𝑖 ~ Normal(𝜇𝑖 , 1) ⚫ 𝜇𝑖 ~ Normal 0, 𝑒 𝛼 2 ⚫ 𝛼 ~ Normal 0,3 ◼ 再パラメータ化したモデル ⚫ 𝑌𝑖 ~ Normal(𝜇𝑖 , 1) ⚫ 𝜇𝑖 = 𝜇𝑖∗ 𝑒 𝛼/2 ⚫ 𝛼 = 3𝛼 ∗ ⚫ 𝜇𝑖∗ ~ Normal 0,1 ⚫ 𝛼 ∗ ~ Normal 0,1 標準正規分布からの サンプリング

5/20 例1）ニールの漏斗 (Neal’s Funnel) ※松浦 (2016, Ch. 9) を参考にした。 ⚫ パラメータの分布 (サンプリング空間) が複雑だと， 事後分布全体から効率よくサンプリングすることができない。 ⚫ 直接推定するパラメータの事後分布はなるべく単純な形にし，それを変形することで 本命のパラメータの事後分布を生成すれば，効率よくサンプリングできる。 𝛼 𝛼 𝜇1 𝜇1 ※松浦 (2016, Ch. 9) の図10.5を改変。

6/20 例2）正規分布 (基本形なので必ず押さえる！) ◼ 再パラメータ化前 ⚫ 𝑌𝑖𝑗 ~ Normal(𝛼𝑖 , 𝜎𝜀 ) ⚫ 𝛼𝑖 ~ Normal 𝜇𝛼 , 𝜎𝛼 ⚫ 𝜇𝛼 ~ Uniform −100,100 ⚫ 𝜎𝛼 ~ Uniform 0,100 ⚫ 𝜎𝜀 ~ Uniform 0,100 ◼ 再パラメータ化後 ⚫ 𝑌𝑖𝑗 ~ Normal(𝜇𝛼 + 𝜎𝛼 𝛼𝑖∗ , 𝜎𝜀 ) ⚫ 𝛼𝑖∗ ~ Normal 0,1 標準正規分布からの サンプリング ⚫ 𝜇𝛼 ~ Uniform −100,100 ⚫ 𝜎𝛼 ~ Uniform 0,100 ⚫ 𝜎𝜀 ~ Uniform 0,100 𝑋 ~ Normal 𝜇, 𝜎 , 𝑍 ~ Normal 0,1 のとき， 𝑋 = 𝜇 + 𝜎𝑍と表せることを利用している。

7/20 例3）対数正規分布 ◼ 再パラメータ化前 ⚫ 𝑌𝑖𝑗 ~ Normal(𝛼𝑖 , 𝜎𝜀 ) ⚫ 𝛼𝑖 ~ LogNormal 𝜇𝛼 , 𝜎𝛼 ⚫ 𝜇𝛼 ~ Uniform −100,100 ⚫ 𝜎𝛼 ~ Uniform 0,100 ⚫ 𝜎𝜀 ~ Uniform 0,100 ◼ 再パラメータ化後 ∗ 𝜇 +𝜎 𝛼 𝛼 𝛼 𝑖 ⚫ 𝑌𝑖𝑗 ~ Normal(𝑒 , 𝜎𝜀 ) 標準正規分布からの ⚫ 𝛼𝑖∗ ~ Normal 0,1 サンプリング ⚫ 𝜇𝛼 ~ Uniform −100,100 ⚫ 𝜎𝛼 ~ Uniform 0,100 ⚫ 𝜎𝜀 ~ Uniform 0,100 𝑋 ~ LogNormal 𝜇, 𝜎 , 𝑊 ~ Normal 𝜇, 𝜎 のとき， 𝑋 = 𝑒 𝑊 と表せることを利用している。

8/20 切断正規分布の 再パラメータ化

9/20 切断正規分布とは ◼ 切断正規分布 (truncated normal distribution) ⚫ 下側または上側の範囲が制限された正規分布。 ⚫ 範囲が(𝑎, 𝑏)の切断正規分布の確率密度関数： 正規分布の確率密度関数 𝑥−𝜇 2 2𝜎2 2𝜋𝜎2 𝑏−𝜇 𝑎−𝜇 Φ −Φ 𝜎 𝜎 1 • 𝑝 𝑥|𝜇, 𝜎, 𝑎, 𝑏 = exp − 𝑎 𝑏 正規分布の(a,b)の範囲の面積 (総和を1に調整する役割) ※Φは標準正規分布の累積分布関数 ⚫ 本発表では簡単のため下限が0の切断正規分布のみを考える。 1 • Normal+ 𝑥|𝜇, 𝜎 = 𝑝 𝑥|𝜇, 𝜎 = 2𝜋𝜎2 exp − Φ 𝑥−𝜇 2 2𝜎2 𝜇 𝜎 0 Normal+に従うパラメータの再パラメータ化を考える

10/20 再パラメータ化の難しさと方針 ◼ 切断正規分布は切断標準正規分布の線形結合では容易に表せない ⚫ 切断の相対位置が𝜇と𝜎に依存するため。 e.g., Normal+ (𝜇, 𝜎)に対応する切断標準正規分布の切断位置は0ではない → 相似な形の切断標準正規分布の表現は複雑 ⚫ ChatGPTに再パラメータ化の方法を聞いたら最初はこの方法を提案されたが， 明らかに間違いであった。 ◼ 累積分布関数の逆関数を使えば一様分布から作れることに気付いた ⚫ ChatGPTにしつこく質問したらこの方法を提案してくれた。 ⚫ (0,1)の一様乱数を累積分布関数の逆関数の引数にいれて乱数生成すると， 使用した関数の元となる確率分布からの乱数を生成できる(結構よく使う)。 e.g., Φ−1 に(0,1)の一様乱数を与えると標準正規分布からの乱数が得られる。 → 切断正規分布の累積分布関数を使って再パラメータ化を試みた。

11/20 最終的に採用したやり方 ◼ 再パラメータ化前 ⚫ 𝑌𝑖𝑗 ~ Normal(𝛼𝑖 , 𝜎𝜀 ) ⚫ 𝛼𝑖 ~ Normal+ 𝜇𝛼 , 𝜎𝛼 ※𝑖は参加者，𝑗は試行 (反復) ※その他の事前分布は簡単のため improperな一様分布を使用。 ◼ 再パラメータ化後 ⚫ 𝑌𝑖𝑗 ~ Normal(𝛼𝑖 , 𝜎𝜀 ) ⚫ 𝛼𝑖 = min 𝜇𝛼 − 𝜎𝛼 ∙ Φ−1 𝛼𝑖∗ ∙ Φ ⚫ 𝛼𝑖∗ ~ Uniform 0,1 𝜇𝛼 𝜎𝛼 , 105 (0,1)の一様分布からのサンプリング

12/20 再パラメータ化前のモデル

13/20 再パラメータ化後のモデル

14/20 解説① 1. (1,0)の一様分布に従う𝛼𝑖∗ を作る。 元の正規分布におけるここの面積 ⚫ 𝛼𝑖∗ ~ Uniform 0,1 ＝1−Φ ⚫ 単純な確率分布からのサンプリング 2. 𝜇𝛼 𝜎𝛼 𝛼𝑖∗ の範囲を 1 − Φ 0−𝜇𝛼 𝜎𝛼 𝜇 = Φ 𝜎𝛼 𝛼 , 1 に制約する 𝜇 ⚫ 𝛼𝑖∗ → 1 − 𝛼𝑖∗ ∙ Φ 𝜎𝛼 𝛼 0 ⚫ 元の正規分布の累積確率に対応させている。 (切断されている部分の累積確率を取らない一様分布を切断) 3. 対応する切断標準正規分布の値を得る。 𝜇 𝜇 𝛼 𝛼 ⚫ 1 − 𝛼𝑖∗ ∙ Φ 𝜎𝛼 → Φ−1 1 − 𝛼𝑖∗ ∙ Φ 𝜎𝛼 𝜇 = −Φ−1 𝛼𝑖∗ ∙ Φ 𝜎𝛼 𝛼 ⚫ Step 2により，切断点より小さな値にはならないことが保証されている。 4. 切断標準正規分布を線形変換してスケールを戻し，式を整理する。 𝜇 ⚫ 𝛼𝑖 = 𝜇𝛼 + 𝜎𝛼 ∙ −Φ−1 𝛼𝑖∗ ∙ Φ 𝜎𝛼 𝛼 𝜇 = 𝜇𝛼 − 𝜎𝛼 ∙ Φ−1 𝛼𝑖∗ ∙ Φ 𝜎𝛼 𝛼

15/20 解説② ◼ このままMCMCすると， 𝛼𝑖∗ がたまに無限大に発散して警告が出る。 (計算上は問題はないがとても鬱陶しい。) ◼ 𝛼𝑖∗ が大きすぎるときには適当な値に置換してやることで警告を防げる。 ⚫ 𝛼𝑖 = min 𝜇𝛼 − 𝜎𝛼 ∙ Φ−1 𝛼𝑖∗ ∙ Φ 𝜇𝛼 𝜎𝛼 , 105

16/20 実際に比較してみよう：Rコード(1/2)

17/20 実際に比較してみよう：Rコード(2/2)

18/20 実際に比較してみよう：結果 再パラメータ化前 再パラメータ化後 • 再パラメータ化前のモデルでは パラメータリカバリも微妙だった。 (𝜇𝛼 のEAPが0.56 [真値は3]) • Rhatはどちらも問題なかった。

19/20 気付いたことなど ◼ Stan上でのPhiの計算がとても不安定 (発散したり誤差が生じたりしやすい) ⚫ なるべくPhiを使わない表現にしたほうがよい。 • 最初は以下のような書き方をしていたが， 𝛼𝑖 が負になるなど計算が不安定だった。 𝛼𝑖 = 𝜇𝛼 − 𝜎𝛼 ∙ Φ−1 Φ − 𝜇𝛼 𝜇𝛼 + 𝛼𝑖∗ 1 − Φ − 𝜎𝛼 𝜎𝛼 ⚫ Phi_approxを代わりに使う手もあるが，近似を許容するかは分析者次第。 (ロジスティック関数を使った近似。詳細はBowling et al., 2009参照。) ◼ 真値によっては再パラメータ化するとかえって収束が悪くなる ⚫ なので必ずしも再パラメータ化したほうが安定するというわけではない。 ⚫ うまく収束しなかったときの1つの解決策として試してみるぐらいがちょうどよさそう。 ◼ target記法で書く場合でも対数尤度の引き算をしなくてよい。 ⚫ The “neglecting the vectorization” errorを未然に防げる。 (Tsukamura & Okada, 2024) ◼ 可逆変換なので周辺尤度の計算には影響しない。 ⚫ bridge samplingするときは自動でヤコビアン調整されるので大丈夫

20/20 引用文献 ⚫ Bowling, S. R., Khasawneh, M. T., Kaewkuekool, S., & Cho, B. R. (2009). A logistic approximation to the cumulative normal distribution. Journal of Industrial Engineering and Management, 2(1), 114–127. https://doi.org/10.3926/jiem.2009.v2n1.p114-127 ⚫ 松浦 健太郎 (2016). StanとRでベイズ統計モデリング 共立出版 ⚫ Tsukamura, Y., & Okada, K. (2024). The “ neglecting the vectorization " error in Stan : erroneous coding practices for computing marginal likelihood and Bayes factors in models with vectorized truncated distributions. Behaviormetrika, 51(2), 635–644. https://doi.org/10.1007/s41237-02400232-7

• https://doi.org/10.3926/jiem.2009.v2n1.p114-127
• https://doi.org/10.1007/s41237-024-00232-7