多変量対数正規分布の計算 (basic statistics of multivariate lognormal distribution)

機械学習・統計

多変量正規分布については色々な解説がありますが、多変量対数正規分布についての結果が知りたいと思ったら意外と(英語でも)解説が見つからなかったので、仕方なく自分で計算しました。メモ代わりにまとめておきます。

Basic statistics(means and (co)variances) of multivariate lognormal distribution are calculated below. Although explanations are in Japanese, you might see mathematical expressions and get information you want. Comments are welcome if you have any questions.

多変量対数正規分布の確率密度関数

多変量対数正規分布の確率密度関数は
P(\bf{y})=\left(\frac{1}{\sqrt{2\pi}}\right)^p \frac{1}{|\bf{\Sigma}|}\frac{1}{y_1y_2\cdots y_p}e^{-\frac{1}{2}(\log \bf{y}-\bf{\mu})^T\bf{\Sigma}^{-1}(\log \bf{y}-\bf{\mu})}
という形になります。ベクトル \bf{\mu}や行列 \bf{\Sigma}は多変量正規分布の平均ベクトル、分散共分散行列に対応しますが、多変量対数正規分布の平均や分散共分散にはなっていないことに注意が必要です。これらのパラメータを用いて、多変量対数正規分布の平均ベクトルと分散共分散行列を求めることが目的になります。

変数変換

1変数の対数正規分布は、 x=(\log y - \mu)/\sigmaと変数変換をすると xについての標準正規分布になります。これに対応する変数変換のを多変量対数正規分布に対しても行うことで、各種統計量の計算を行います。
今、 \bf{\Sigma}の固有値 \lambda_kと固有ベクトル \bf{u}_kが分かったとします。このとき \bf{\Sigma}は
\Sigma_{ij} = \sum_k \lambda_k u_{ki}u_{kj}
および \bf{\Sigma}^{-1}は
\Sigma^{-1}_{ij} = \sum_k \frac{1}{\lambda_k} u_{ki}u_{kj}
と書けます（証明略）。ここで、 u_{ki}はk番目の固有ベクトルの第i成分を表します。これを利用して、変数変換を
t_i = \sum_j \frac{1}{\sqrt{\lambda_i}}u_{ij}(\log y_j - \mu_j)
とします。多変量の積分の変数変換にはヤコビアンを計算する必要があり、
\frac{d(t_1, \cdots, t_p)}{d(y_1, \cdots, y_p)} = \left| \begin{matrix} \frac{dt_1}{dy_1} & \cdots & \frac{dt_1}{dy_p} \\ \vdots & & \vdots \\ \frac{dt_p}{dy_1} & \cdots & \frac{dt_p}{dy_p} \\ \end{matrix} \right| = \left| \begin{matrix} \frac{u_{11}}{\sqrt{\lambda_1}y_1} & \cdots & \frac{u_{1p}}{\sqrt{\lambda_1}y_p} \\ \vdots & & \vdots \\ \frac{u_{p1}}{\sqrt{\lambda_p}y_1} & \cdots & \frac{u_{pp}}{\sqrt{\lambda_p}y_p} \\ \end{matrix} \right|
となります。行列式の性質を使って計算を進めると
\frac{d(t_1, \cdots, t_p)}{d(y_1, \cdots, y_p)} = \frac{1}{y_1\cdots y_p}\frac{1}{\sqrt{\lambda_1 \cdots \lambda_p}}\left|U|
となります( U_{ij} = u_{ij})。Uは直交行列なので、行列式は1となり、結局ヤコビアンは
\frac{d(t_1, \cdots, t_p)}{d(y_1, \cdots, y_p)} = \frac{1}{y_1\cdots y_p}\frac{1}{\sqrt{\lambda_1 \cdots \lambda_p}}
となります。yからtへ変換する際には逆変換となり、その場合、ヤコビアンは逆数となります。

平均

y_iの平均値は
E(y_i)=\int dy_1 \cdots dy_p \left(\frac{1}{\sqrt{2\pi}}\right)^p \frac{1}{|\bf{\Sigma}|}\frac{y_i}{y_1\cdots y_p}e^{-\frac{1}{2}(\log \bf{y}-\bf{\mu})^T\bf{\Sigma}^{-1}(\log \bf{y}-\bf{\mu})}
となります。先ほどの変数変換を行うと、確率密度関数の部分は無相関な標準正規分布となり、 y_iは
y_i = \e\left({\mu_i + \sum_j u_{ji}\sqrt{\lambda_j}t_j}\right)
となります。ここで Uは直交行列であり、逆行列が U^{-1} = U^Tであるという性質を使っています。従って、
E(y_i)=\int dt_1 \cdots dt_p \left(\frac{1}{\sqrt{2\pi}}\right)^p \e\left({\mu_i + \sum_j u_{ji}\sqrt{\lambda_j}t_j}\right) e^{-\frac{1}{2}\bf{t}^2
となります。積分を計算するためには各 t_jについて平方完成すればよく、eの肩に乗っている項のうち関連するもののみを見ると、
-\frac{1}{2}t_j^2 + u_{ji}\sqrt{\lambda_j}t_j = -\frac{1}{2}(t_j - u_{ji}\sqrt{\lambda_j})^2 + \frac{1}{2}u_{ji}^2 \lambda_j
とすることができます。 -\frac{1}{2}(t_j - u_{ji}\sqrt{\lambda_j})^2の項は t_jで積分すると規格化定数 \frac{1}{\sqrt{2\pi}}とキャンセルして消えるので、結局残るのは
E(y_i)= \e \left(\mu_i + \frac{1}{2} \sum_j u_{ji}^2 \lambda_j \right)
になります。ここで、
\Sigma_{ij} = \sum_k \lambda_k u_{ki}u_{kj}
という式を思い出すと、
\frac{1}{2} \sum_j u_{ji}^2 \lambda_j = \frac{1}{2} \Sigma_{ii} = \frac{1}{2} \sigma_i^2
であることが分かります。結局、多変量対数正規分布の平均値は
E(y_i)= \e \left(\mu_i +\frac{1}{2} \sigma_i^2)
と計算できることが分かりました。

分散共分散

平均の次は分散(共分散)を計算します。
共分散の計算はやや面倒ですが、考え方は平均のときと同じものを使うことができます。
Cov(y_i,y_j)=\int dy_1 \cdots dy_p \left(\frac{1}{\sqrt{2\pi}}\right)^p \frac{1}{|\bf{\Sigma}|}\frac{(y_i-\bar{y}_i)(y_j-\bar{y}_j)}{y_1\cdots y_p}e^{-\frac{1}{2}(\log \bf{y}-\bf{\mu})^T\bf{\Sigma}^{-1}(\log \bf{y}-\bf{\mu})}
平均値 \bar{y}_i, \bar{y}_jは既に求まっているのでそれを代入し、平均の計算と同様に標準正規分布への変数変換を行うと
Cov(y_i,y_j)= \int dt_1 \cdots dt_p \left(\frac{1}{\sqrt{2\pi}}\right)^p \left( e^{\mu_i + u_{ki}\sqrt{\lambda_k}t_k} - e^{\mu_i + \frac{1}{2}\sigma_i^2} \right) \left( e^{\mu_j + u_{kj}\sqrt{\lambda_k}t_k} - e^{\mu_j + \frac{1}{2}\sigma_j^2} \right \big) e^{-\frac{1}{2}\bf{t}^2}
ただし、ここでスペース節約のため、アインシュタインの縮約記法（単一の項内で重複する添字については和を取る）を用いています。すなわち、eの肩に載っている項のkについては和を取るという意味になります。
これらの項のうち、積分後に生き残る（ゼロにならない）のは(eの肩に)tの1乗が含まれる項です。また生き残る項についての積分の結果は、平均の計算のときと同様に、元々の確率密度関数から来る t^2項と合わせた平方完成による「おつり」の項です。この計算を地道に行うと
Cov(y_i,y_j) = e^{\mu_i + \mu_j} \left( e^{\frac{1}{2}\lambda_k u_{ki}^2 + \frac{1}{2}\lambda_k u_{kj}^2 + \lambda_k u_{ki}u_{kj}} - e^{\frac{1}{2}(\sigma_i^2 + \sigma_j^2)}\right)
となり、平均値のときと同様、 \lambda_k u_{ki}^2 = \sigma_i^2、 \lambda_k u_{ki}u_{kj} = \sigma_{ij}となるので、結局
Cov(y_i,y_j) = e^{\mu_i + \frac{1}{2}\sigma_i^2 + \mu_j + \frac{1}{2}\sigma_j^2}\left( e^{\sigma_{ij}} - 1 \right) = E(y_i)E(y_j) \big( e^{\sigma_{ij}} - 1\big)
となります。