確率統計
多変量正規分布サンプリング; Multi-variate Normal Distribution Sampling

概要

---

多変量正規分布サンプリング
分散共分散行列\(\boldsymbol{\Sigma} \), 平均\(\boldsymbol{\mu} \)の\(N\)次元正規分布(多変量正規分布)は以下の式で表現される確率分布。

$$ f(\boldsymbol{x}) = \frac{1}{\sqrt{| 2 \pi \boldsymbol{\Sigma} | }} exp \left( -\frac{1}{2} (\boldsymbol{x} - \boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} (\boldsymbol{x} - \boldsymbol{\mu}) \right) $$


多変量正規分布に従う確率変数の生成は以下の手順で行う。
1.分散共分散行列から下三角行列\(L=(l_{ij})_{1\leq i, j \leq N}\)を生成
2.標準正規分布に従う確率変数\(Z_1, Z_2, \ldots , Z_N\)を生成
3.以下の式に基づいて確率変数\( \boldsymbol{x} \)を求める
\( \begin{eqnarray} x_1 &=& \mu_1 + l_{11} Z_1 \\ x_2 &=& \mu_2 + l_{21} Z_1 + l_{22} Z_2 \\ \vdots \\ x_N &=& \mu_N + l_{N1} Z_1 + l_{N2} Z_2 + \ldots + l_{NN} Z_N \\ \end{eqnarray} \)

ソースコード

---

namespace ExRandom.MultiVariate {
    public class MultiVariateNormalRandom : Random<double> {
        readonly Continuous.NormalRandom nd;
        readonly int dim;
        readonly double[] mu_vector;
        readonly double[][] lower_tri_matrix;

        public MultiVariateNormalRandom(MT19937 mt, double[,] cov_matrix, double[] mu_vector = null){
            if(mt == null) {
                throw new ArgumentNullException();
            }

            if(cov_matrix == null) {
                throw new ArgumentException();
            }

            if(mu_vector != null && mu_vector.Length != cov_matrix.GetLength(0)){
                throw new ArgumentException();
            }

            double[][] l;

            CholeskyDecomp(cov_matrix, out l);

            for(int i = 0, j; i < dim; i++) {
                for(j = 0; j <= i; j++) {
                    if(Double.IsNaN(l[i][j])) { 
                        throw new ArgumentException();
                    }
                }
            }

            this.nd = new Continuous.NormalRandom(mt);
            this.dim = cov_matrix.GetLength(0);
            this.lower_tri_matrix = l;

            this.mu_vector = (mu_vector != null) ? mu_vector : new double[dim];
        }

        public override Vector<double> Next() {
            double[] z = nd.Next(dim), x = new double[dim];

            for(int i = 0, j; i < dim; i++) { 
                x[i] = mu_vector[i];
                
                for(j = 0; j <= i; j++) {
                    x[i] += lower_tri_matrix[i][j] * z[j];
                }
            }

            return new Vector<double>(x);
        }

        private static void CholeskyDecomp(double[,] m, out double[][] l) {
            if(m.GetLength(0) <= 0 || m.GetLength(0) != m.GetLength(1)) {
                throw new ArgumentException();
            }

            int i, j, k, n = m.GetLength(0);

            for(i = 1; i < n; i++) {
                for(j = 0; j < i; j++) {
                    if(m[i, j] != m[j, i]) { 
                        throw new ArgumentException();
                    }
                }
            }

            double ld, lld;
            double[] d = new double[n];
            
            l = new double[n][];

            for(i = 0; i < n; i++) { 
                l[i] = new double[i + 1];
            }

            l[0][0] = 1;
            d[0] = m[0, 0];

            for(i = 0; i < n; i++) {
                for(j = 0; j < i; j++) {
                    lld = m[i, j];
                    for(k = 0; k < j; k++) {
                        lld -= l[i][k] * l[j][k] * d[k];
                    }
                    l[i][j] = lld / d[j];
                }

                ld = m[i, i];
                for(k = 0; k < i; k++) { 
                    ld -= l[i][k] * l[j][k] * d[k];
                }
                l[i][j] = 1;
                d[i] = ld;
            }

            for(j = 0; j < n; j++) { 
                d[j] = Math.Sqrt(d[j]);
            }

            for(i = 0; i < n; i++) {
                for(j = 0; j <= i; j++) {
                    l[i][j] *= d[j];
                }
            }
        }
    }
}

関連項目

---

メルセンヌ・ツイスタ
各種確率分布サンプリング基本クラス