1

我正在尝试编写一个 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

任何帮助将不胜感激！

4

1 回答 1

0

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

问题是将 iid N(0,1) 随机数的向量转换为目标多元法线。检查教科书，您需要协方差矩阵 S 的 Cholesky 分解 A

S = AA'

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

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

于 2013-05-14T15:46:16.647 回答