はじめに

多次元ガウス分布では、相関の変化を検知して異常検知を行うスパースモデリングが存在する。そこで、計数の分布に対しても同じようなモデリングが出来るのではないかと考え、多次元ベルヌーイ分布を調べてみたが、思ったよりも相関構造が複雑で、難しいことがわかった。

参考

読書日記: 読了：Teugels (1990) 多変量ベルヌーイ分布をどうやって表現するか
Teugels, J.L. (1990) Some representation of the multivariate bernoulli and binomial distributions. Journal of Multivariate Analysis. 32(2), 256-268.

多次元ベルヌーイ分布

$\{ X_i \}$ をベルヌーイ分布に従う変数列とする。周辺確率分布は次で与えられるとする。

\begin{aligned}
p_i &= \mathrm{Pr}(X_i=1) = E[X_i], \\
q_i &= \mathrm{Pr}(X_i=0) = 1 - p_i = 1 - E[X_i]
\end{aligned}

次に多次元ベルヌーイ分布を考える。多次元ベルヌーイ分布を表現するのに、次の表記を導入する。

p_{k_1, k_2, \cdots, k_n} = \mathrm{Pr}(X_1=k_1, X_2=k_2, \cdots, X_n=k_n), \ k_i\in \{ 0, 1\}

導入として、2次元のときから考えてみる。2次元のベルヌーイ分布を考える上で自然なパラメータは次の3つである。

\begin{aligned}
p_1 &= \mathrm{Pr}(X_1=1) \\
p_2 &= \mathrm{Pr}(X_2=1) \\
\sigma_{12} &= E[(X_1 - p_1)(X_2 - p_2)] = E[X_1 X_2] - p_1 p_2
\end{aligned}

これらのパラメータを用いて、同時分布$p_{k_1k_2}$を表現する。

\begin{aligned}
p_{11} &= \mathrm{Pr}(X_1=1, X_2=1) \\
&= \sum_{k_1, k_2}k_1k_2\mathrm{Pr}(X_1=k_1, X_2=k_2) \\
&= E[X_1X_2] \\
&= p_1p_2+\sigma_{12} \\
&:= \mu_{12} \\ 
p_{10} &= \mathrm{Pr}(X_1=1, X_2=0) \\
&= \sum_{k_1, k_2}k_1(1-k_2)\mathrm{Pr}(X_1=k_1, X_2=k_2) \\
&= E[X_1(1-X_2)] \\
&= p_1 - \mu_{12} \\
&= p_1q_2 - \sigma_{12} \\
p_{01} &= \mathrm{Pr}(X_1=0, X_2=1) \\
&= \sum_{k_1, k_2}(1-k_1)k_2\mathrm{Pr}(X_1=k_1, X_2=k_2) \\
&= E[(1-X_1)X_2] \\
&= p_2 - \mu_{12} \\
&= q_1p_2 - \sigma_{12} \\
p_{00} &= \mathrm{Pr}(X_1=0, X_2=0) \\
&= \sum_{k_1, k_2}(1-k_1)(1-k_2)\mathrm{Pr}(X_1=k_1, X_2=k_2) \\
&= E[(1-X_1)(1-X_2)] \\
&= 1- p_1 - p_2 + \mu_{12} \\
&= q_1q_2 + \sigma_{12} \\
\end{aligned}

当たり前だが、$\sum p_{k_1k_2} = p_{00}+p_{10}+p_{01}+p_{11}=1$が成り立つ。

同時分布を縦に並べて、ベクトルとしてみる。

\begin{aligned}
\boldsymbol{p}^{(2)}
:= \left(
    \begin{array}{ccc}
      p_{00} \\
      p_{10} \\
      p_{01} \\
      p_{11} 
    \end{array}
  \right)
&= \left(
    \begin{array}{ccc}
      q_1 q_2 & -q_2 & -q_1 & 1 \\
      p_1 q_2 & q_2 & -p_1 & -1 \\
      q_1 p_2 & -p_2 & q_1 & -1 \\
      p_1 p_2 & p_2 & p_1 & 1
    \end{array}
  \right)
\left(
    \begin{array}{ccc}
      1 \\
      0 \\
      0 \\
      \sigma_{12} 
    \end{array}
  \right)  \\
&= \left(
    \begin{array}{ccc}
      1 & -1 & -1 & 1 \\
      0 & 1 & 0 & -1 \\
      0 & 0 & 1 & -1 \\
      0 & 0 & 0 & 1
    \end{array}
  \right)
\left(
    \begin{array}{ccc}
      1 \\
      p_1 \\
      p_2 \\
      \mu_{12} 
    \end{array}
  \right)
\end{aligned}

これらの表現はテンソル積でスマートに書ける。

\begin{aligned}
\boldsymbol{p}^{(2)}
&= \left(
    \begin{array}{ccc}
      q_2 & -1 \\
      p_2 & 1
    \end{array}
  \right)
 \otimes \left(
    \begin{array}{ccc}
      q_1 & -1 \\
      p_1 & 1
    \end{array}
  \right)
\left(
    \begin{array}{ccc}
      1 \\
      0 \\
      0 \\
      \sigma_{12} 
    \end{array}
  \right) \\
& =\left(
    \begin{array}{ccc}
      1 & -1 \\
      0 & 1
    \end{array}
  \right)
 \otimes \left(
    \begin{array}{ccc}
      1 & -1 \\
      0 & 1
    \end{array}
  \right)
\left(
    \begin{array}{ccc}
      1 \\
      p_1 \\
      p_2 \\
      \mu_{12} 
    \end{array}
  \right)
\end{aligned}

3次元の場合は次のように書ける。

\begin{aligned}
\boldsymbol{p}^{(2)}
&= \left(
    \begin{array}{ccc}
      q_3 & -1 \\
      p_3 & 1
    \end{array}
  \right)
 \otimes \left(
    \begin{array}{ccc}
      q_2 & -1 \\
      p_2 & 1
    \end{array}
  \right)
 \otimes \left(
    \begin{array}{ccc}
      q_1 & -1 \\
      p_1 & 1
    \end{array}
  \right)
\left(
    \begin{array}{ccc}
      1 \\
      0 \\
      0 \\
      \sigma_{12} \\
      0 \\
      \sigma_{13} \\
      \sigma_{23} \\
      \theta \\
    \end{array}
  \right)
\end{aligned}

ここで、 $\sigma_{ij} = E[(X_i - p_i)(X_j - p_j)], \theta=E[(X_1 - p_1)(X_2 - p_2)(X_3 - p_3)]$ である。

これを$n$次元に拡張する。多次元ベルヌーイ分布の実現値$(k_1, k_2, \cdots, k_n)$を表現するために次の量を導入する。

k = 1 + \sum_{i=1}^n k_i 2^{i-1}, \ k_i\in \{ 0, 1\}, \ 1\leq k \leq 2^n

$k$は、$(k_1, k_2, \cdots, k_n)$と1対1対応する。確率も次のように書く。

p^{(n)}_k = p_{k_1, k_2, \cdots, k_n}

多次元ベルヌーイ分布の確率質量関数は次のように書ける。

p^{(n)}_k = \mathrm{Pr}(X_1=k_1, X_2=k_2, \cdots, X_n=k_n)
= E\left[\prod_{i=1}^n X_i^{k_i} \bar X_i^{1-k_i} \right]

ここで$\bar X_i=1-X_i$である。

2行1列ベクトルのテンソル積に成り立つ公式、

\left( 
\left(
    \begin{array}{ccc}
      a_n \\
      b_n 
    \end{array}
\right)\otimes
\left(
    \begin{array}{ccc}
      a_{n-1} \\
      b_{n-1} 
    \end{array}
\right)\otimes
\cdots\otimes
\left(
    \begin{array}{ccc}
      a_{1} \\
      b_{1} 
    \end{array}
\right)
 \right)_k
= \prod_{i=1}^k  a_i^{1-k_i} b_i^{k_i}

を用いると、確率ベクトルは次のように書ける。

\boldsymbol{p}^{(n)} = E\left[ 
\left(
    \begin{array}{ccc}
      \bar X_n \\
      X_n 
    \end{array}
\right)\otimes
\left(
    \begin{array}{ccc}
      \bar X_{n-1} \\
      X_{n-1} 
    \end{array}
\right)\otimes
\cdots\otimes
\left(
    \begin{array}{ccc}
      \bar X_{1} \\
      X_{1} 
    \end{array}
\right)
 \right]

さらに、これをモーメントを使用して書き直す。モーメントベクトルのk成分は次のように定義する。

\begin{aligned}
\mu_k^{(n)} = E \left[ \prod_{i=1}^n X_i^{k_i} \right]
= E\left[\left( 
\left(
    \begin{array}{ccc}
      1 \\
      X_n 
    \end{array}
\right)\otimes
\left(
    \begin{array}{ccc}
      1 \\
      X_{n-1} 
    \end{array}
\right)\otimes
\cdots\otimes
\left(
    \begin{array}{ccc}
      1 \\
      X_{1} 
    \end{array}
\right)
 \right)_k
\right]
\end{aligned}

また、中心モーメントも定義する。

\begin{aligned}
\sigma_k^{(n)} = E \left[ \prod_{i=1}^n (X_i-p_i)^{k_i} \right]
= E\left[\left( 
\left(
    \begin{array}{ccc}
      1 \\
      Y_n 
    \end{array}
\right)\otimes
\left(
    \begin{array}{ccc}
      1 \\
      Y_{n-1} 
    \end{array}
\right)\otimes
\cdots\otimes
\left(
    \begin{array}{ccc}
      1 \\
      Y_{1} 
    \end{array}
\right)
 \right)_k
\right]
\end{aligned}

ここで、$Y_i = X_i-p_i$である。

論文では次の定理が示されている。

\boldsymbol{p}^{(n)} = 
\left(
    \begin{array}{ccc}
      1 & -1 \\
      0 & 1
    \end{array}
\right)^{\otimes n} \boldsymbol{\mu}^{(n)} ,\\
\boldsymbol{\mu}^{(n)} = 
\left(
    \begin{array}{ccc}
      1 & 1 \\
      0 & 1
    \end{array}
\right)^{\otimes n} \boldsymbol{p}^{(n)} ,\\
\boldsymbol{p}^{(n)} = 
\left(
    \begin{array}{ccc}
      q_n & -1 \\
      p_n & 1
    \end{array}
\right)\otimes
\left(
    \begin{array}{ccc}
      q_{n-1} & -1 \\
      p_{n-1} & 1
    \end{array}
\right)
 \otimes \cdots \otimes\left(
    \begin{array}{ccc}
      q_1 & -1 \\
      p_1 & 1
    \end{array}
  \right)

 \boldsymbol{\sigma}^{(n)} ,\\
\boldsymbol{\sigma}^{(n)} = 
\left(
    \begin{array}{ccc}
      1 & 1 \\
      -p_n & q_n
    \end{array}
\right)\otimes
\left(
    \begin{array}{ccc}
      1 & 1 \\
      -p_{n-1} & q_{n-1}
    \end{array}
\right)
 \otimes \cdots \otimes\left(
    \begin{array}{ccc}
      1 & 1 \\
      -p_1 & q_1
    \end{array}
  \right)
\boldsymbol{p}^{(n)} ,\\

$\mu_1^{(n)}=1$だから、一つ目の式から$\boldsymbol{p}^{(n)}$の自由度は$2^n-1$となる。三つ目の式からの自由度勘定も考えてみる。$\boldsymbol{\sigma}^{(n)}$のうち、$\sigma_1^{(n)}=1$で、任意のiで$k_i=1, k_{i以外}=0$のようなときは$\sigma_k^{(n)}=E[Y]=0$となる。したがって$n+1$の自由度が引かれ、$2^n-n-1$となる。これに$(p_1,\cdots,p_n)$のn個の自由度があるため、結局$2^n-1$となる。

$n$次元ガウス分布の場合、分布は平均 $\boldsymbol{\mu}$、分散共分散行列 $\Sigma$ で特徴づけられる。したがって自由度は$n+n(n+1)/2$となる。これは高次のモーメントも$\boldsymbol{\mu}$と分散共分散行列 $\Sigma$で表されることを示している。特に2変数間の偏相関は$\Sigma$の非対角成分で表される。一方で、多次元ベルヌーイ分布は高次のモーメントも独立な自由度を持っているため、$2^n-1$の自由度を持つ。ガウス分布は自由度にかなり強い制限がかかっていると考えられる。

ベルヌーイ分布の偏相関を考えてみる。

\begin{aligned}
\mathrm{Pr}(X_1=k_1, X_2=k_2 | X_3=k_3, \cdots, X_n=k_n)
&= \frac{\mathrm{Pr}(X_1=k_1, X_2=k_2, \cdots, X_n=k_n)}{\mathrm{Pr}(X_3=k_3, \cdots, X_n=k_n)} \\
&\propto E\left[X_1^{k_1} \bar X_1^{1-k_1} X_2^{k_2} \bar X_2^{1-k_2} \prod_{i=3}^n X_i^{k_i} \bar X_i^{1-k_i} \right]
\end{aligned}

これから、$X_1$と$X_2$の直接相関がないための条件は次のようになる。

\begin{aligned}
X_1 \perp X_2 | X_3,\cdots,X_n \Leftrightarrow X_1 \perp \mathrm{otherｓ} \lor X_2 \perp \mathrm{otherｓ}
\end{aligned}

他の変数を通じて関連しあっているため、ガウス分布のように条件が$\Sigma^{-1}_{12}=0$とはならない。相関の推定も、全変数に対する相関を調べる必要がある。

多次元二項分布

結論だけ。

ベルヌーイ分布$\mathrm{Pr}(X_1=k_1, X_2=k_2, \cdots, X_n=k_n)$から$m$個のサンプルを取得したとする。$X_i$の$m$個のサンプルの和を$S_{im}=X_{i1}+\cdots+X_{im}$とする。$\{S_{im}\},i=1,\cdots,n$の同時分布を次のように書く。

\begin{aligned}
\pi^{(m)}_{r_1,r_2,\cdots,r_n} = \mathrm{Pr}(S_{1m}=r_1, S_{2m}=r_2, \cdots, S_{nm}=r_n)
\end{aligned}

$(r_1,r_2,\cdots,r_n)$と1対1対応する量$k$を次のように与える。

k = 1 + \sum_{i=1}^n r_i(m+1)^{i-1}, \ r_i\in \{0,1,\cdots,m\}, 1\leq k \leq (m+1)^n

ベクトル表記を次のように導入する。

\begin{aligned}
_n\boldsymbol{q}^{(m)} &:= (_n q_1^{(m)}, _n q_2^{(m)}, \cdots, _n q_{(m+1)^n}^{(m)}) \\
_n q_k^{(m)} &:= \pi^{(m)}_{r_1,r_2,\cdots,r_n}, \ 1\leq k \leq (m+1)^n
\end{aligned}

$_n\boldsymbol{q}^{(m)}$を$p_i$と$\sigma_k^{(n)}$で表すことが出来る。

\begin{aligned}
_n\boldsymbol{q}^{(m)} &= V_m^{(n)}V_{m-1}^{(n)}\cdots V_{1}^{(n)} \\
V_i^{(n)} &= \sum_{k=1}^{2^n}\sigma_k^{(n)}R_{n,i}^{(k)}\otimes \cdots \otimes R_{1,i}^{(k)} \\
R_{m,i}^{(k)} &= 
  \begin{cases}
    Q_{i,m} & (k_i=0)\\
    J_m  & (k_i=1)
  \end{cases}\\
Q_{i,m}&=\left(
    \begin{array}{cccccc}
      q_i & 0 & 0 & \cdots & 0 & 0 \\
      p_i & q_i & 0 & \cdots & 0 & 0 \\
      0 & p_i & q_i & \cdots & 0 & 0 \\
      \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\
      0 & 0 & 0 & \cdots & p_i & q_i \\
      0 & 0 & 0 & \cdots & 0 & p_i
    \end{array}
\right)\\
J_m &= \left(
    \begin{array}{cccccc}
      -1 & 0 & 0 & \cdots & 0 & 0 \\
      1 & -1 & 0 & \cdots & 0 & 0 \\
      0 & 1 & -1 & \cdots & 0 & 0 \\
      \cdots & \cdots & \cdots & \cdots & \cdots & \cdots \\
      0 & 0 & 0 & \cdots & 1 & -1 \\
      0 & 0 & 0 & \cdots & 0 & 1
    \end{array}
\right)
\end{aligned}

ここで、$Q_{i,m}$および$J_m$は$(m+1)\times m$行列。