从 fortran 中的条件多元正态分布中提取答案

【问题标题】：Draw from a conditional multivariate normal distribution in fortran从 fortran 中的条件多元正态分布中提取
【发布时间】：2013-05-08 03:59:24
【问题描述】：

我正在尝试编写一个 fortran 子例程，以根据另一个子空间的状态从多元正态分布中提取子样本。基本上：

(x1, x2)' ~ N( (mu1, mu2)', 西格玛)

其中协方差矩阵Sigma可以划分为四个子矩阵

西格玛=（S11，S12；S21，S22）

教科书和维基百科告诉我 x1 在 x2=a 上的条件分布是：

x1|x1=a ~ N(mu, Sigma*)

在哪里

mu = mu1 + S12 * S22^-1 * (a - mu2)

西格玛* = S11 - S12 * S22^-1 * S21

当在 R 中编写它时，它就像一个魅力。在 Fortran 中没有那么多。

SUBROUTINE dCondMVnorm ( DIdx, NDraw, Sigma, NSigma, Mu, TMCurr)

IMPLICIT NONE

INTEGER            :: I, NSigma, NDraw, INFO
INTEGER            :: DIdx(NDraw), NIdx(NSigma-NDraw), AllIdx(NSigma)
LOGICAL            :: IdxMask(NSigma)
DOUBLE PRECISION   :: Sigma11(NDraw, NDraw), Sigma22(NSigma-NDraw,NSigma-NDraw)
DOUBLE PRECISION   :: Sigma(NSigma,NSigma)
DOUBLE PRECISION   :: Sigma12minv22(NSigma-NDraw,NDraw), iSigma22(NSigma-NDraw,NSigma-NDraw)
DOUBLE PRECISION   :: RandNums(NDraw), Dummy1(NDraw), MeanDiff(NSigma-NDraw )
DOUBLE PRECISION   :: TMcurr(NSigma), Mu(NSigma)

! create the indeces to _not_ draw from (NIdx)
IdxMask = .FALSE.
IdxMask(DIdx) = .TRUE.
AllIdx = (/ (I, I=1, NSigma) /)
NIdx = pack( AllIdx, .NOT. IdxMask)

Sigma11 = Sigma( DIdx, DIdx)
Sigma22 = Sigma( NIdx, NIdx)
iSigma22 =0.0D0
DO I = 1, NSigma-NDraw
  iSigma22(I,I) = 1.0D0
END DO
CALL DPOSV( 'U', NSigma-NDraw,NSigma-NDraw, Sigma22, NSigma-NDraw, iSigma22, NSigma-NDraw, INFO)
CALL DGEMM ( 'N', 'N', NDraw, NSigma-NDraw, NSigma-NDraw, 1.0D0, Sigma(DIdx,NIdx), NDraw, &
   iSigma22, NSigma-NDraw, 0.0D0, Sigma12minv22, NDraw )

CALL DGEMM ( 'N', 'N', NDraw, NDraw, NSigma-NDraw, -1.0D0, Sigma12minv22, NDraw, &
   Sigma(NIdx,DIdx), NSigma-NDraw, +1.0D0, Sigma11, NDraw)
CALL DPOTRF( 'U', NDraw, Sigma11, NDraw, INFO)
DO I = 1, NDraw-1
  Sigma11(I+1:NDraw,I) = 0.0D0
END DO
! now Sigma11 actually is the cholesky decomposition of the matrix Sigma*
MeanDiff = TMcurr(NIdx) - Mu(NIdx)
CALL DGEMV( 'N', NDraw, NSigma-NDraw, 1.0D0, Sigma12minv22, NDraw, MeanDiff, 1, 0.0D0, Dummy1(1), 1)

! sorry, this one is self written and returns NDraw random numbers that are i.i.d. N(0,1) using Marsaglia's algorithm
CALL getzig(RandNums, NDraw)
CALL DGEMV( 'N', NDraw, NDraw, 1.0D0, Sigma11, NDraw, RandNums(1), 1, 1.0D0, Dummy1(1), 1)

TMcurr(DIdx) = Dummy1
END SUBROUTINE dCondMVnorm

所以我现在构建这个（它是我正在处理的更大模块的一部分）从 R 中调用它

CovMat <- diag(4)
CovMat[1:3,2:4] <- CovMat[1:3,2:4] + diag(3)*.5
CovMat[2:4,1:3] <- CovMat[2:4,1:3] + diag(3)*.5
CovMat[3:4,1:2] <- CovMat[3:4,1:2] + diag(2)*.2
CovMat[1:2,3:4] <- CovMat[1:2,3:4] + diag(2)*.2
library(MASS)
dyn.load("TM_Updater.so")
testMat2 <- matrix(NA,0,4)
for (a in seq(500) ){
  y <- mvrnorm(1,rep(0,2), CovMat[3:4,3:4])
  x <- .Fortran("dCondMVnorm", as.integer(c(1,2)),as.integer(2), CovMat, as.integer(4), c(0.0,0.0,0.0,0.0), c(0.0,0.0,y))[[6]]
  testMat2 <- rbind(testMat2, c(x[1:2],y) )
}
dyn.unload("TM_Updater.so")
cov(testMat2)

这会返回

> cov(testMat2)
        [,1]      [,2]      [,3]        [,4]
[1,] 1.179618573 0.4183372 0.1978489 0.002156081
[2,] 0.418337156 0.8317497 0.4891746 0.204091537
[3,] 0.197848928 0.4891746 0.9649001 0.498660858
[4,] 0.002156081 0.2040915 0.4986609 1.032272666

显然，[1,1] 的协方差太高了，无论我多久运行一次（或运行多长时间）都是如此。我错过了什么？ Fortran 计算的协方差矩阵与手工计算的协方差矩阵相匹配，平均值也是如此……精度不同的一些问题？

另外，DGEMV 有一个奇怪的地方，你需要给出向量 y 的确切起始地址（参见对 DGEMV 的最后一次调用）（如纪录片中所称）才能得到

y := alpha A *x + beta * y, beta != 0

任何帮助将不胜感激！

【问题讨论】：

标签： r fortran distribution

【解决方案1】：

我感到很尴尬，但为了将来参考，这个错误应该让所有人都能看到。

问题是转换 i.i.d 的向量。目标多元正态分布的 N(0,1) 个随机数。查课本需要协方差矩阵S的cholesky分解A

S = AA'

请注意，我们感兴趣的是下三角矩阵，而不是我计算的上三角矩阵。

解决方案：在对 DGEMV 的最后一次调用中，将“N”更改为“T”或在对 DPOSV 的调用中计算“L”三角形，并在以下几行中重写上三角形的归零。

【讨论】：