多変量正規分布にしたがう乱数の生成


多変量正規分布にしたがう乱数を生成したい場合のお話です。正規分布は条件付き確率も正規分布なので Gibbs サンプラーを容易に構築できるのですが,互いに独立な乱数が生成できるのであればそれにこしたことはありません。ということで互いに独立な多変量正規乱数を生成できるようにします。

を各次元の要素が独立に標準単変量正規分布 にしたがう 次元の乱数ベクトルであるとします。言い換えれば は多変量正規分布 ( は零ベクトル, は単位行列) にしたがっています。これから となる変換をほどこしたときに が目的分布 ( は共分散行列) にしたがうように を定めます。そうすれば は明らかに にしたがいます。

による変換により,次のようにして の平均が で共分散行列が となることがわかります。

したがって となるように を定めれば良いことになります。ここで は実対称行列なので対角化可能で となるような直行行列 および固有値を対角成分に持つ行列 が存在します。さらに は共分散行列なので,半正定値行列です。つまり固有値はすべて 0 以上の実数です。 とすれば とすることができます。そして は対角行列なので となります。したがって となり, とすれば良いことがわかります。

さて実際にこのような を構築するのには,たとえばヤコビ法やコレスキー分解を用います。いずれも適用できる条件がありますが,共分散行列はいずれの条件も満たしています。ここでは気分的にヤコビ法を採用します[ とする時点でコレスキー分解で良いのですが,最初その発想にいたらなくて対称行列からヤコビ法しか知らなかったのでこれを選んだというのは内緒。">A]

let sqrtsumsq x y =
   if abs x > abs y then
      let r = y / x
      in abs x * sqrt (1.0 + r * r)
   elif y <> 0.0 then
      let r = x / y
      in abs y * sqrt (1.0 + r * r)
   else
      0.0

let diag n = Array2D.init n n (fun i j -> if i = j then 1.0 else 0.0)

let epsilon = pown 2.0 (-52)
let jacobi matrix =
   let n = Array2D.length1 matrix
   let m = n - 1
   let eigenvalues = Array2D.copy matrix
   let eigenvectors = diag n
   let findMax () =
      seq {
         for i = 0 to m do
            for j = 0 to m do
               if i <> j then
                  yield (i, j), abs eigenvalues.[i, j]
      }
      |> Seq.maxBy snd
   let loop = ref true
   while !loop do
      let (p, q), max = findMax ()
      if max < epsilon then
         loop := false
      else
         let app = eigenvalues.[p, p]
         let aqq = eigenvalues.[q, q]
         let apq = eigenvalues.[p, q]
         let t = 0.5 * (app - aqq)
         let ss = 0.5 * (1.0 - abs t / sqrtsumsq apq t)  // sin^2
         let cc = 1.0 - ss  // cos^2
         let s = if apq * t > 0.0 then -sqrt ss else sqrt ss  // sin
         let c = sqrt cc  // cos
         let sc = s * c  // sin * cos
         for i = 0 to m do
            let api = eigenvalues.[p, i]
            let aqi = eigenvalues.[q, i]
            eigenvalues.[p, i] <- api * c - aqi * s
            eigenvalues.[q, i] <- api * s + aqi * c
         for i = 0 to m do
            eigenvalues.[i, p] <- eigenvalues.[p, i]
            eigenvalues.[i, q] <- eigenvalues.[q, i]
         eigenvalues.[p, p] <- app * cc - 2.0 * apq * sc + aqq * ss
         eigenvalues.[q, q] <- aqq * cc + 2.0 * apq * sc + app * ss
         eigenvalues.[p, q] <- 0.0
         eigenvalues.[q, p] <- 0.0
         for i = 0 to m do
            let aip = eigenvectors.[i, p]
            let aiq = eigenvectors.[i, q]
            eigenvectors.[i, p] <- aip * c - aiq * s
            eigenvectors.[i, q] <- aip * s + aiq * c
   eigenvalues, eigenvectors

事前に必要となる標準正規乱数の生成は,たとえば Box-Muller 変換を用います。

let random = System.Random ()
let rec boxMuller () = seq {
   let r = sqrt <| -2.0 * log (random.NextDouble ())
   let theta = 2.0 * System.Math.PI * (random.NextDouble ())
   yield r * cos theta
   yield r * sin theta
   yield! boxMuller ()
}

これらを用いて多変量正規分布にしたがう乱数を生成します。

let multiply a b =
   let m = Array2D.length1 a
   let k = Array2D.length2 a
   let n = Array2D.length2 b
   let p = Array2D.zeroCreate m n
   for i = 0 to m - 1 do
      for j = 0 to n - 1 do
         let mutable sum = 0.0
         for t = 0 to k - 1 do
            sum <- sum + a.[i, t] * b.[t, j]
         p.[i, j] <- sum
   p

let multiplyVector a x =
   let m = Array2D.length1 a
   let n = Array.length x
   let p = Array.zeroCreate m
   for i = 0 to m - 1 do
      let mutable sum = 0.0
      for t = 0 to n - 1 do
         sum <- sum +  a.[i, t] * x.[t]
      p.[i] <- sum
   p

let multinormal mu sigma =
   let n = Array.length mu
   let eigenvalues, eigenvectors = jacobi sigma
   let d = Array2D.map sqrt eigenvalues
   let q = multiply eigenvectors d
   let z = boxMuller () |> Seq.take n |> Seq.toArray
   Array.map2 (+) mu (multiplyVector q z)

ということでいつでも多変量正規分布にしたがう互いに独立な乱数が生成できるようになりました。

脚注

  1. 本当は正定値行列を とする時点でコレスキー分解で良いのですが,最初その発想にいたらなくて対称行列からヤコビ法しか知らなかったのでこれを選んだというのは内緒。 []