他のサブスペースの状態を条件とする多変量正規分布からサブサンプルを描画するための 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
どんな助けでも大歓迎です!