Draw from a conditional multivariate normal distribution in fortran

Question

I am trying to write a fortran subroutine to draw a subsample from a multivariate normal distribution conditional on the state of the other subspace. Basically:

(x1, x2)' ~ N( (mu1, mu2)', Sigma)

where the covariance matrix Sigma can be partitioned in the four submatrices

Sigma=( S11, S12; S21, S22)

Textbooks & Wikipedia tell me that the conditional distribution of x1 on x2=a is:

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

where

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

Sigma* = S11 - S12 * S22^-1 * S21

When writing this up in R it works like a charm. In Fortran not so much.

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

So I now build this (it is part of a larger module I am working on) call this from R using

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)

and this returns

> 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

clearly, the covariance of [1,1] is much too high and it is that way no matter how often (or for how long) I run it. What am I missing? The covariance matrix calculated by Fortran matches the one calculated by hand, as do the means... some issues with different accuracies?

Plus there's this weirdness with the DGEMV that you need to give the exact starting address (see last call to DGEMV) of the vector y (as it is called in the documentary) in order to get

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

Any help would greatly be appreciated!

Joe W Joe W · Accepted Answer · 2013-05-14T15:46:16

I feel embarrassed, but for future reference this blunder shall remain available for all to see.

The problem is converting a vector of i.i.d. N(0,1) random numbers to the target multivariate normal. Checking the textbooks you need the cholesky decomposition A of the covariance matrix S

S = AA'

Note that it is the lower triangular matrix we are interested in, not the upper that I calculated.

Solution: in the last call to DGEMV change 'N' to 'T' or calculate the 'L' triangle in the call to DPOSV and rewrite the zeroing out of the upper triangle in the following lines.

Draw from a conditional multivariate normal distribution in fortran

1 Answers