问具有Math.net的多元正态分布
EN

Stack Overflow用户

提问于 2012-06-30 11:47:50

回答 1查看 1.7K关注 0票数 0

我插入Math.net做一些抽样，我只是不明白它的工作方式.我在这里错过了什么？

let mm  =  Double.DenseMatrix.Identity(2)
let ida =  Double.DenseMatrix.Identity(2)
let idb =  Double.DenseMatrix.Identity(2)
let generator = MatrixNormal(mm, ida, idb)
generator.Density(mm)

我得到了

System.ArgumentOutOfRangeException: Matrix dimensions must agree.
Parameter name: x
   at MathNet.Numerics.Distributions.MatrixNormal.Density(Matrix`1 x) in c:\TeamCity\buildAgent\work\392bcd0e1411b00f\src\Numerics\Distributions\Multivariate\MatrixNormal.cs:line 241
   at <StartupCode$FSI_0079>.$FSI_0079.main@()
Stopped due to error

奇怪的是，当我查看github源代码时，就会抛出

  public double Density(Matrix<double> x)        { 
       if (x.RowCount != _m.RowCount || x.ColumnCount != _m.ColumnCount) 
       {                
           throw Matrix.DimensionsDontMatch<ArgumentOutOfRangeException>(x, _m, "x");                         
       }

编辑

在重新启动一切之后，它对2*2矩阵的均值起作用，但对2*1 (调整列方差后)不起作用。非常奇怪的是，在定义上有一些维度检查。然而，错误消息是在调用时出现的。可能是检查错误行方差和列方差，并在调用时使用正确的。所有这些都突出了富类型检查的强大优势。

对于好奇的人来说，这里是一个多元的实现。不过还没查过。

let sampleNormal = 
   let rnd = new MersenneTwister()
   fun () ->
      let rec randomNormal () = 
          let u1, u2 = rnd.NextDouble(),rnd.NextDouble()
          let r = sqrt (-2. * (log u1))
          let theta = 2. * System.Math.PI * u2  
          seq { yield r * sin theta
                yield r * cos theta 
                yield! randomNormal() }
      randomNormal ()


let generate covar = 
   let chol = Double.Factorization.DenseCholesky(covar)
   let a = chol.Factor
   fun () -> let v = vector ( sampleNormal() |> Seq.take(covar.ColumnCount) |> List.ofSeq )
             a * v
//generate covar  

let generatecovar = generate covar
let generaten n covar = Seq.init n (fun _ -> generatecovar ())

编辑

Cholesky对于完全相关的输入失败，这是可以的。

let mapply (m:Generic.Matrix<float>) f = m.IndexedEnumerator() |> Seq.iter(fun (i,j,v) -> m.[i,j] <- f v ); m

let generate (covar:Generic.Matrix<float>) = 
   let R = if covar.Determinant() = 0. then   // we want covar = R.RT  // C = U D1/2.D1/2 U' = (RT.QT) Q.R = RT.RTT  
               let u, d, vt = let t = covar.Svd(true) in t.U(), t.W(), t.VT()  
               let A = (mapply d sqrt) * u.Transpose()                    
               let qr = A.QR() in qr.R.Transpose()                                 
            else
               let chol = covar.Cholesky()
               chol.Factor
   fun () -> let v = vector ( sampleNormal() |> Seq.take(covar.ColumnCount) |> List.ofSeq )
             R * v