1

他のサブスペースの状態を条件とする多変量正規分布からサブサンプルを描画するための Fortran サブルーチンを作成しようとしています。基本的：

(x1, x2)' ~ N( (mu1, mu2)', シグマ)

ここで、共分散行列シグマは 4 つの部分行列に分割できます

シグマ=( S11, S12; S21, S22)

教科書とウィキペディアによると、x2=a での x1 の条件付き分布は次のようになります。

x1|x1=a ~ N( μ, シグマ*)

どこ

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 によって計算された共分散行列は、手段と同様に、手で計算されたものと一致します...異なる精度の問題はありますか?

さらに、取得するためにベクトル y の正確な開始アドレス (DGEMV への最後の呼び出しを参照) を指定する必要がある DGEMV には、この奇妙さがあります (ドキュメンタリーで呼び出されているように)。

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

どんな助けでも大歓迎です！

4

1 に答える 1

0

恥ずかしい思いをしますが、今後の参考のために、この失敗はすべての人が見ることができるように残します.

問題は、iid N(0,1) 乱数のベクトルをターゲットの多変量法線に変換することです。教科書を確認すると、共分散行列 S のコレスキー分解 A が必要です

S = AA'

計算したのは上三角行列ではなく、下三角行列であることに注意してください。

解決策: DGEMV への最後の呼び出しで、'N' を 'T' に変更するか、DPOSV への呼び出しで 'L' 三角形を計算し、次の行で上の三角形のゼロ設定を書き直します。

于 2013-05-14T15:46:16.647 に答える