ひっさびさのブログ記事更新です。
最近技術系の記事を書いてなかったなと思って、お正月休みを利用して書いてみようと思います。
今回紹介するのは、Stanで総和が0になるように制約する方法です。
え、そんなのに工夫がいるのか？と思われるかもですが、実はいろいろややこしい話があったりするのでそれをまとめてみました。あと、cmdstan2.36ではsum_to_zero_vectorという新しい型のベクトルが登場しました。その使い方についても解説します。
◆総和を0に制約する
N次元ベクトルxで総和を0に制約したいとき、N-1次元のベクトルyを自由パラメータとして、x[N]の値を-sum(y)とする、すなわちx = append_row(y, -sum(y))とする方法が使われることが多かったと思います。しかし、この方法では、N-1次元まではyのままで、x[N]だけが新たに構成されてしまうことから、非対称な事前分布からサンプリングしていることになります。たとえば総和を0に制約した回帰効果（分散分析などがそれ）を見たいときに、パラメータの事前分布が非対称になってしまうのは、推定上問題が生じることがあります。
具体的にどのように非対称になるのかをシミュレーションで確認してみましょう。
|
1 2 3 4 5 6 7 8 |
n <- 5 m <- 100000 y <- rnorm(m*(n-1),0,1) |> array(dim=c(m,n-1)) x <- array(NA,dim=c(m,n)) for(i in 1:m){ x[i,] <- c(y[i,],-sum(y[i,])) } |
上のように、まず、独立に正規分布からN-1次元のベクトルを10万個生成します。そのあと、x[N]を-sum(y)で計算した場合の分布をみてみましょう。
|
1 2 |
x |> apply(2,sd) x |> cor() |
|
1 2 3 4 5 6 7 8 9 |
> x |> apply(2,sd) [1] 0.9993724 0.9975755 0.9982909 0.9979529 1.9987791 > x |> cor() [,1] [,2] [,3] [,4] [,5] [1,] 1.000000000 0.007637069 -0.0008764760 -0.0030917309 -0.5018216 [2,] 0.007637069 1.000000000 -0.0047945315 0.0049451104 -0.5029853 [3,] -0.000876476 -0.004794531 1.0000000000 0.0005583468 -0.4968980 [4,] -0.003091731 0.004945110 0.0005583468 1.0000000000 -0.5004823 [5,] -0.501821644 -0.502985269 -0.4968979909 -0.5004823457 1.0000000 |
このように、N-1次元目までは標準偏差が1、相関は0となりますが、N次元目だけ標準偏差が2、ほかのベクトルの要素との相関は―0.5と非常に高い負の相関が出てしまいます。
本来、総和が0であるという制約のみをあたえて対称な事前分布からxを生成する場合、全体的に負の相関が生じますが、上の方法では1つのベクトルの要素にだけ強い負の相関が生じ、ほかが無相関という偏った事前分布から生成されることになってしまいます。このことから、従来の方法では適切なベイズ推定にならない場合があります。
これを回避するため、yをそのまま活用するのではなく、全体的に変換を行う必要があります。この方法は複数考えられますが、Stanのマニュアルではこのページに書かれている方法などが提案されています（なお、最新のStanマニュアルに書かれている方法ではうまく計算ができない…と思うので、上のサイトにある変換方法を紹介します）。
R上で、sum_to_zero()関数を以下のように定義します。
|
1 2 3 4 5 6 7 |
sum_to_zero <- function(y){ n <- length(y)+1 norm2 <- sum(y)/(sqrt(n)+n) fill_val = norm2 - sum(y)/sqrt(n) x <- c(y,fill_val)-norm2 return(x) } |
この関数を使って、ベクトルの総和が0になるように変換してみましょう。
|
1 2 3 4 |
y <- c(1,2,3,4) x <- sum_to_zero(y) x sum(x) |
|
1 2 3 4 |
> x [1] -0.381966 0.618034 1.618034 2.618034 -4.472136 > sum(x) [1] 0 |
たしかに、総和が0になっています。
それでは、この変換方法では事前分布がどのようになるのかを確認してみます。
|
1 2 3 4 5 6 7 8 9 10 |
n <- 5 m <- 100000 y <- rnorm(m*(n-1),0,1) |> array(dim=c(m,n-1)) x <- array(NA,dim=c(m,n)) for(i in 1:m){ x[i,] <- y[i,] |> sum_to_zero() } x |> apply(2,sd) |> mean() x |> cor() |
|
1 2 3 4 5 6 7 8 9 |
> x |> apply(2,sd) |> mean() [1] 0.894951 > x |> cor() [,1] [,2] [,3] [,4] [,5] [1,] 1.0000000 -0.2474729 -0.2478553 -0.2548746 -0.2486143 [2,] -0.2474729 1.0000000 -0.2496913 -0.2464573 -0.2562136 [3,] -0.2478553 -0.2496913 1.0000000 -0.2521535 -0.2483653 [4,] -0.2548746 -0.2464573 -0.2521535 1.0000000 -0.2482977 [5,] -0.2486143 -0.2562136 -0.2483653 -0.2482977 1.0000000 |
このように、すべてのベクトルの要素の標準偏差と相関が一定になります。理論値的には、sum_to_zero変換したベクトルの各要素の標準偏差は、sqrt((n-1(/n)となります。n=5の場合は、理論値は0.8944となり、シミュレーション結果はそれとほぼ一致しています。
◆Stanのsum_to_zero_vector型
cmdstan2.36から、sum_to_zero_vectorという新しい型が登場しました。このベクトルは、上で説明したように対称な事前分布になるようにうまいこと変換してくれるので、特に何も考えずに総和が0になるようなベクトルを簡単に推定することができます。
ただ、事前分布を仮定するときに1つ注意が必要なのが、sum_to_zeroベクトルはベクトルの次元数に応じて小さくなってしまう（上の理論値を参照）点です。よって、事前分布は次元数に応じてやや大きめ(sqrt(N/(N-1)))に設定しておく必要があります。そこでStanでsum_to_zero_vectorを使うとき、事前分布には次のように設定します。
|
1 |
x ~ normal(mu, sigma*sqrt(N/(N-1))); |
ただ、推定効率上、実際は次のようにしたほうがよいです（数値的には同じ）。
|
1 |
x ~ normal(mu, sigma*inv(sqrt(1-inv(N)))); |
今後のバージョンでは、zero_sum_normal_lpdfという新しい分布が実装される予定のようです。これを使えば、標準偏差も自動的に次元数に応じて調整してくれます。
◆simplex型をsoftmax変換で作る場合も使える
総和が1に制約されたベクトルはsimplexといいます。Stanではsimplex型がすでにありますが、モデリングの都合によってはN-1次元の独立した自由パラメータから、x = softmax(append_row(0, y))というように変換したいこともあるでしょう。実はsimplexに変換するときもこの方法では非対称な事前分布になってしまいます。
確認してみましょう。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
softmax <- function(x){ exp(x)/sum(exp(x)) } n <- 5 m <- 100000 y <- rnorm(m*(n-1),0,1) |> array(dim=c(m,n-1)) x <- array(NA,dim=c(m,n)) for(i in 1:m){ x[i,] <- c(0,y[i,]) |> softmax() } x |> apply(2,sd) x |> cor() |
|
1 2 3 4 5 6 7 8 9 |
> x |> apply(2,sd) [1] 0.07507438 0.16307381 0.16317828 0.16362190 0.16277274 > x |> cor() [,1] [,2] [,3] [,4] [,5] [1,] 1.0000000 -0.1197285 -0.1114080 -0.1158373 -0.1131450 [2,] -0.1197285 1.0000000 -0.3156867 -0.3161687 -0.3123369 [3,] -0.1114080 -0.3156867 1.0000000 -0.3174824 -0.3156983 [4,] -0.1158373 -0.3161687 -0.3174824 1.0000000 -0.3167631 [5,] -0.1131450 -0.3123369 -0.3156983 -0.3167631 1.0000000 |
やはり、最初の要素だけが標準偏差と相関が、他の要素と違っています。
この問題も、上で紹介したsum_to_zero関数を使えば解決します。
|
1 2 3 4 5 6 7 8 9 10 |
n <- 5 m <- 100000 y <- rnorm(m*(n-1),0,1) |> array(dim=c(m,n-1)) x <- array(NA,dim=c(m,n)) for(i in 1:m){ x[i,] <- y[i,] |> sum_to_zero() |> softmax() } x |> apply(2,sd) x |> cor() |
|
1 2 3 4 5 6 7 8 9 |
> x |> apply(2,sd) [1] 0.1621395 0.1630671 0.1629523 0.1624981 0.1630252 > x |> cor() [,1] [,2] [,3] [,4] [,5] [1,] 1.0000000 -0.2520866 -0.2497746 -0.2464691 -0.2470806 [2,] -0.2520866 1.0000000 -0.2489635 -0.2494968 -0.2519974 [3,] -0.2497746 -0.2489635 1.0000000 -0.2509913 -0.2519283 [4,] -0.2464691 -0.2494968 -0.2509913 1.0000000 -0.2511963 [5,] -0.2470806 -0.2519974 -0.2519283 -0.2511963 1.0000000 |
もちろん、Stan内蔵のsimplex型を使えばもともと問題はないのですが、自分で何とかしたい場合などはこの方法がいいでしょう。