Рисовать из условного многомерного нормального распределения в Фортране

Question

Я пытаюсь написать подпрограмму на фортране, чтобы нарисовать подвыборку из многомерного нормального распределения, зависящего от состояния другого подпространства. В принципе:

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

где ковариационная матрица Sigma может быть разделена на четыре подматрицы

Сигма =( S11, S12; S21, S22)

Учебники и Википедия говорят мне, что условное распределение x1 по x2=a:

x1 | x1 = a ~ N (мю, сигма *)

где

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

Сигма * = S11 - S12 * S22^-1 * S21

При написании этого в R это работает как шарм. В Фортране не так много.

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] слишком велика, и это происходит независимо от того, как часто (или как долго) я ее запускаю. Что мне не хватает? Матрица ковариации, вычисленная Фортраном, совпадает с матрицей, вычисленной вручную, как и средства... некоторые проблемы с различной точностью?

Плюс в DGEMV есть та странность, что вам нужно дать точный начальный адрес (см. Последний вызов DGEMV) вектора y (как его называют в документальном фильме), чтобы получить

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

Любая помощь будет принята с благодарностью!

Answer 1

Я чувствую смущение, но для дальнейшего использования этот промах будет доступен для всеобщего обозрения.

Задача состоит в преобразовании вектора случайных чисел iid N(0,1) в целочисленную многомерную нормаль. При проверке учебников нужно разложить холески A ковариационной матрицы S

S = AA '

Обратите внимание, что интересует нижняя треугольная матрица, а не верхняя, которую я рассчитал.

Решение: в последнем вызове DGEMV замените 'N' на 'T' или вычислите треугольник 'L' в вызове DPOSV и перепишите обнуление верхнего треугольника в следующих строках.