« 軽井沢駅前の喫茶店の店員さんに渡されたパズルを解く (ためのRパッケージを作った) | メイン | 読了：Ahmad & Khan (2019) 量質混在データのクラスタリング手法レビュー »
2019年11月28日 (木)
しばらく前に、$n$個の二値変数の同時分布を少数のパラメータでうまく表現しなければならない、という用事があり、さんざっぱら悩む羽目になった。地べたを這いずり回るような仕事をしているもので、時々そういう変な用事が出現するのである。なんかさあ！AIがどうとかさあ！ディープラーニングがどうとかさあ！そういう華やかな仕事をしたいもんだよね！ (嘘です。仕事せず寝ているのが理想だ)
この話題、専門の方には簡単なんだろうけど、私にはかなり難しくて、途中でかなりの混乱に陥った。恥ずかしながら当初は、多変量正規分布と同じパラメータ数があれば十分だろう、つまり、周辺割合が$n$個、2変数間の連関を表すパラメータが$n(n-1)/2$個あれば十分だろう、と思い込んでいたのである。なんとなく多次元IRTやカテゴリカル因子分析モデルが頭にあって、個々の二値変数の背後に正規潜在変数があるというのを自明視していたような気がする。いやー、おかげでずいぶん回り道をした。
冷静に考えると、必要なパラメータ数は$2^n-1$個である。なぜ多変量正規分布よりもパラメータ数が増える？？とまたまた混乱したのだが、お世話になっている先生に意外な角度からのアドバイスを頂いたりして(さすがプロの研究者はちがう...)、数日かけて頭を冷やし、素人ながらもなんとなく腑に落ちてきた。普段あまり意識していなかったけど、多変量正規性というのは結構きつい制約なんですね。反省。
というわけで、個人的には学ぶところ多かったが、これを考える事情をもたらした案件のほうは中止になってしまい、私の苦闘は丸ごと無駄に終わった。どうもそういう人生みたいだ。
Teugels, J.L. (1990) Some representation of the multivariate bernoulli and binomial distributions. Journal of Multivariate Analysis. 32(2), 256-268.
上記の都合により、勉強のつもりで読んだ論文。
多変量ベルヌーイ分布、つまり複数のベルヌーイ変数の分布について考える。
まずは2変数の場合について考えよう。以下、$p_{00}=P(X_1 = 0, X_2 = 0)$という風に書くことにする。
2変数の場合のパラメータは3つである。その3つは、$p_{00}, p_{10}, p_{01}, p_{11}$のうちどれか3つだと考えてもいいし(当然ながら4つの和は1だから)、ふたつの周辺分布$p_1 = E(X_1) = P(X_1 = 1), p_2$と、$\sigma_{12} = E(X_1-p_1)(X_2-p_2)$の計3個だ、と考えてもいい。
[ちょっとわかりにくいけど、$p_{00}$の添字は値を示し、$p_1$や$\sigma_{12}$の添字は変数番号を示している]
$q_1 = 1-p_1$という風に書くとして、
$p_{00} = q_1 q_2 + \sigma_{12}$
$p_{10} = p_1 q_2 - \sigma_{12}$
$p_{01} = q_1 p_2 - \sigma_{12}$
$p_{11} = p_1 p_2 + \sigma_{12}$
と書ける。行列だと、$\mathbf{p}^{(2)} = [p_{00}, p_{10}, p_{01}, p_{11}]^T$として、
$\displaystyle \mathbf{p}^{(2)} = \left[\begin{array}{cc} q_2 & -1 \\ p_2 & 1 \end{array} \right] \otimes \left[\begin{array}{cc} q_1 & -1 \\ p_1 & 1 \end{array} \right] \left[\begin{array}{c} 1 \\ 0 \\ 0 \\ \sigma_{12} \end{array}\right]$
と書ける。
[ええええ？ という感じだけど... $\otimes$はクロネッカー積を表していて、この例だと、右側の行列に$q_2$を掛けたのを左上、$-1$を掛けたのを右上、$p_2$を掛けたのを左下、$1$を掛けたのを右下においた4x4の行列ができる。その右から掛けている縦ベクトルは2行目と3行目が0だから、この行列の2列目と3列目は無視するとして、1列目は上から$q_1 q_2, p_1 q_2, q_1 p_2, p_1 p_2$になるし、4列目は上から$1, -1, -1, 1$になるから、なるほど、行列にする前と比べてつじつまが合っている]
あるいは、$\sigma_{12}$のかわりに$\mu_{12} = E X_1 X_2 = \sigma_{12}+p_1 p_2 = p_{11}$を使ってもよい。このとき、
$p_{00} = 1-p_1-p_2+\mu_{12}$
$p_{10} = p_1 -\mu_{12}$
$p_{01} = p_2 - \mu_{12}$
$p_{11} = \mu_{12}$
と書ける。[こっちのほうがなんとなく親しみがわきますね]
行列だと
$\displaystyle \mathbf{p}^{(2)} = \left[\begin{array}{cc} 1 & -1 \\ 0 & 1 \end{array} \right] \otimes \left[\begin{array}{cc} 1 & -1 \\ 0 & 1 \end{array} \right] \left[\begin{array}{c} 1 \\ p_1 \\ p_2 \\ \mu_{12} \end{array}\right]$
となる。[ああ... クロネッカー積というのを持ち出してきた理由がやっとわかった]
3変数に拡張してみよう。
新たに$\theta = E(X_1 - p_1)(X_2 - p_2)(X_3 - p_3)$を導入する。パラメータは$p_1,p_2,p_3,\sigma_{12},\sigma_{13},\sigma_{23},\theta$の7つになり、行列で書くと
$\mathbf{p}^{(3)} = \left[\begin{array}{cc} q_3 & -1 \\ p_3 & 1 \end{array} \right] \otimes \left[\begin{array}{cc} q_2 & -1 \\ p_2 & 1 \end{array} \right] \otimes \left[\begin{array}{cc} q_1 & -1 \\ p_1 & 1 \end{array} \right] \left[ \begin{array}{c} 1 \\ 0 \\ 0 \\ \sigma_{12} \\ 0 \\ \sigma_{13} \\ \sigma_{23} \\ \theta \end{array} \right]$
となる。
[なんだか騙されたような気分だ... なんでこういう風になるのか、いまいち腑に落ちないんだけど、とにかく先に進もう]
$\mathbf{p}^{(n)}$へと一般化しよう。
$\mathbf{p}^{(n)}$は長さ$2^n$のベクトルである。いま、$k_i \in \{0,1\}$として、$(k_1, k_2, \ldots, k_n)$と一対一に対応するbinary expansion形式の変数
$k = 1 + \sum_{i=1}^{n} k_i 2^{i-1}$
を考えよう。$p_{k_1,k_2,\ldots,k_n}$ は$p_k^{(n)}$と書きかえられる。
[なにをほざいておるのかというと、$(k_1, k_2, \ldots, k_n)$に順に$1, 2, 4, 8, \ldots$という重みを振った加重和に1を足して$k$と呼びましょうということである。要するに、$\mathbf{p}^{(n)}$の要素番号を$k$としましょうってことね]
以下では$Y_i = 1-X_i$と書く。[原文では$\bar{X}_i$なんだけど、わかりにくいので勝手に書き換えた]
以下のように書ける。
$p_k^{(n)} = P(\cap_{i=1}^n [X_i = k_i]) = E(\prod_{i=1}^n X_i^{k_i} Y_i^{1-k_i})$
さて、一般に$2 \times 1$ベクトルのクロネッカー積は次の性質を持つ：
$\left[ \left[ \begin{array}{c} a_n \\ b_n \end{array} \right] \otimes \left[ \begin{array}{c} a_{n-1} \\ b_{n-1} \end{array} \right] \otimes \cdots \otimes \left[ \begin{array}{c} a_1 \\ b_1 \end{array} \right] \right]_{k} = \prod_{i=1}^n a_i^{1-k_i} b_i^{k_i}, \ \ 1 < k < 2^n$
[なにをぬかしておるのかというと、たとえば$n=2$のとき、
$\left[ \begin{array}{c} a_2 \\ b_2 \end{array} \right] \otimes \left[ \begin{array}{c} a_1 \\ b_1 \end{array} \right] = \left[ \begin{array}{c} a_2 a_1 \\ a_2 b_1 \\ b_2 a_1 \\ b_2 b_2 \end{array} \right]$
だよね？ で、たとえば$k=2$(つまり$k_1=1, k_2=0$)番目の要素をみると、これは$a_1^{1-k_1} b_1^{k_1} \times a_2^{1-k_2} b_2^{k_2}=a_1^0 b_1^1 \times a_2^1 b_2^0$だよね？ つまり、クロネッカー積の$k$番目の値がほしかったら、その$k$を$(k_1, k_2, \ldots, k_n)$に戻し、$i=1, \ldots, n$について、$k_i$が0だったら$a_i$, $1$だったら$b_i$を拾ってきて掛ければいいんだよ、ということであろう。そういわれればそうなんだけど、よくもまあ、こんな難しい書き方を...]
$a_i = Y_i, b_i = X_i$とすると、
$\mathbf{p}^{(n)} = E \left( \left[ \begin{array}{c} Y_n \\ X_n \end{array} \right] \otimes \left[ \begin{array}{c} Y_{n-1} \\ X_{n-1} \end{array} \right] \otimes \cdots \otimes \left[ \begin{array}{c} Y_1 \\ X_1 \end{array} \right] \right)$
となる。
さて。$\mathbf{\mu}^{(n)}$, $\mathbf{\sigma}^{n}$を次のように定義する。それぞれの$k$番目の要素について、
$\mu_k^{(n)} = E \left( \prod_{i=1}^{n} X_i^{k_i} \right) = E \left( \left[ \begin{array}{c} 1 \\ X_n \end{array} \right] \otimes \left[ \begin{array}{c} 1 \\ X_{n-1} \end{array} \right] \otimes \cdots \otimes \left[ \begin{array}{c} 1 \\ X_1 \end{array} \right] \right)_k$
$\sigma_k^{(n)} = E \left( \prod_{i=1}^{n} (X_i-p_i)^{k_i} \right) = E \left( \left[ \begin{array}{c} 1 \\ Y_n \end{array} \right] \otimes \left[ \begin{array}{c} 1 \\ Y_{n-1} \end{array} \right] \otimes \cdots \otimes \left[ \begin{array}{c} 1 \\ Y_1 \end{array} \right] \right)_k$
[ぴんとこないので$n=2, k=4$(つまり$k_1 = 1, k_2=1$)について考えると、
$\mu_4^{(2)} = E( X_1^{k_1} X_2^{k_2} ) = E( X_1 X_2 ) $
つまりこれは最初の練習に出てきた$\mu_{12}$である。
$\sigma_4^{(2)} = E ( (X_1-p_1)^{k_1} (X_2-p_2)^{k_2} ) = E (X_1-p_1)(X_2-p_2)$
つまりこれは最初の練習に出てきた$\sigma_{12}$である。はいはい]
以下の定理を得ることができます。[証明は省略する。すいません、私ただのしがない労働者ですもんで]
定理1-1.
$\mathbf{p}^{(n)} = \left[ \begin{array}{cc} 1 & -1 \\ 0 & 1 \end{array} \right]^{\otimes n} \mathbf{\mu}^{(n)}$
かつ
$\mathbf{\mu}^{(n)} = \left[ \begin{array}{cc} 1 & 1 \\ 0 & 1 \end{array} \right]^{\otimes n} \mathbf{p}^{(n)}$
パラメータ数は$2^n-1$個となる。$\mathbf{\mu}^{(n)}$の長さは$2^n$だが、$\mu_1^{(n)}=1$だから。
定理1-2.
$\mathbf{p}^{(n)} = \left[ \begin{array}{cc} q_n & -1 \\ p_n & 1 \end{array} \right] \otimes \left[ \begin{array}{cc} q_{n-1} & -1 \\ p_{n-1} & 1 \end{array} \right] \otimes \cdots \otimes \left[ \begin{array}{cc} q_1 & -1 \\ p_1 & 1 \end{array} \right] \mathbf{\sigma}^{(n)}$
かつ
$\mathbf{\sigma}^{(n)} = \left[ \begin{array}{cc} 1 & 1 \\ -p_n & q_n \end{array} \right] \otimes \left[ \begin{array}{cc} 1 & 1 \\ -p_{n-1} & q_{n-1} \end{array} \right] \otimes \cdots \otimes \left[ \begin{array}{cc} 1 & 1 \\ -p_1 & q_1 \end{array} \right] \mathbf{p}^{(n)}$
パラメータ数はどうなるか。$\mathbf{\sigma}^{(n)}$の長さは$2^n$だが、$k_1 + k_2 + \cdots + k_n = 1$の箇所は0になるから、$2^n -n - 1$個。これが依存性を表しているわけだ。これに加えて$p_i$が$n$個あるから、結局パラメータ数は$2^n-1$である。[あー！なるほどねえ！]
[ここからは、多変量ベルヌーイ分布の確率母関数の導出。メモは省略]
なお、分割表の分析には対数線型モデルが使われることも多い。たとえば$n=3$の飽和モデルは、クロネッカー積を使えば次のように書ける。$v_{ijk} = \log p_{ijk}$のベクトルを$\mathbf{v}^{(3)}$として、
$\mathbf{v}^{(3)} = \left[ \begin{array}{cc} 1 & 1 \\ 1 & -1 \end{array} \right]^{\otimes 3} \mathbf{\lambda}^{(3)}$
ただし$\mathbf{\lambda}^{(3)} = (\mu, \lambda_1, \lambda_2, \lambda_{12}, \lambda_3, \lambda_{13}, 0, \phi)^T$。
[ほんまかいな？ というわけで、展開してみると...
$\left[\begin{array}{c} v_{000} \\ v_{100} \\ v_{010} \\ v_{110} \\ v_{001} \\ v_{101} \\ v_{011} \\ v_{111} \end{array} \right] = \left[ \begin{array}{cccccccc} +1 & +1 & +1 & +1 & +1 & +1 & +1 & +1 \\ +1 & -1 & +1 & -1 & +1 & -1 & +1 & -1 \\ +1 & +1 & -1 & -1 & +1 & +1 & -1 & -1 \\ +1 & -1 & -1 & +1 & +1 & -1 & -1 & +1 \\ +1 & +1 & +1 & +1 & -1 & -1 & -1 & -1 \\ +1 & -1 & +1 & -1 & -1 & +1 & -1 & +1 \\ +1 & +1 & -1 & -1 & -1 & -1 & +1 & +1 \\ +1 & -1 & -1 & +1 & -1 & +1 & +1 & -1 \end{array} \right] \left[ \begin{array}{c} \mu \\ \lambda_1 \\ \lambda_2 \\ \lambda_{12} \\ \lambda_3 \\ \lambda_{13} \\ 0 \\ \phi \end{array} \right]$
右辺左側の正方行列は、1列目から順に、全平均、要因$i$の主効果、$j$の主効果, $i \times j$の交互作用効果、$k$の効果、$i \times k$の交互作用効果、$j \times k$の交互作用効果、$i \times j \times k$の二次交互作用効果になっている(というか、L8直交表の左に1の列を付け加えた形になっている)]
これを最初で示した例と見比べると、パラメータ数はどちらも7だけど、解釈が異なる。
我々の$\mathbf{p}^{(n)} = A_n \mathbf{\sigma}^{(n)}$という書き方だと、$\theta=0$とすることなく$\sigma_{12} = \sigma_{13} = \sigma_{23} = 0$とおける。つまり、対数線型モデルのように階層的構造を考える必要がない。
さらに、我々の定式化だと、ひとつ以上の$p_{ijk}$が0であってもなんの問題もない。対数線型モデルでは結構難しくなる。
云々...
... ここまでが論文の前半。後半は多変量二項分布の話になり、いま関心ないのでパス。
習わぬお経を無理矢理読むようなつもりでメモをとったが、それなりに勉強になった。クロネッカー積を使うとこんなに綺麗に書ける、っていうところに感心した。世の中には頭のいい奴がいるものだ。
論文：データ解析(2018-) - 読了：Teugels (1990) 多変量ベルヌーイ分布をどうやって表現するか