Question: Help! Multivarite numerical integration issues with evalf(Int)

Hi.  I'm trying to compute a posterior distribution (Bayesian inference), which requires integration.  I have set up the integral and have attempted to solve it numerically using evalf(Int...) with Maple 12 Worksheet interface and the command-line version.  When executing with Maple 12, I eventually (hours later) get a error - "connection to kernel lost".  When executing with the command-line, eventually the process seems to terminate.  The information about memory used stops updating and nothing happens.  Here is my code:

To keep things focused, i'm omitting the part of the code where I populate the vectors mu and sigma, and the parameters sigm[G], AlphaA, AlphaB, BetaA, BetaB. These computations run fine.  Each element of these vectors is a floating point number. There is also a  matrix of  floating point correlation coefficients called rho that is  already populated.

with(Statistics):
with(LinearAlgebra):

MeanVector := Vector(8, [mu[Field], mu[Stream], mu[Field], mu[Stream], mu[F3], mu[S3], mu[R3], mu[B3]]):
SDVector := Vector(8, [sigma[Field], sigma[Stream], sigma[Field], sigma[Stream], sigma[F3], sigma[S3], sigma[R3], sigma[B3]]):
VarVector := SDVector.Transpose(SDVector):
Sigma := Matrix(8, 8):
X := Vector(8, symbol = x):
for i to 8 do for j to 8 do Sigma[i, j] := rho[i, j]*VarVector[i, j] end do: end do:
a1:= X-MeanVector:
a2:=Transpose(a1):
a3:=1/Sigma:
a4:= exp(-1/2*a2.a3.a1):
MNormDist := 1/((2*Pi)^4*sqrt(Determinant(Sigma)))*a4:

Dist := y^(-1+AlphaA)*(1-y)^(BetaA-1)/Beta(AlphaA, BetaA):
Dist1 := z^(-1+AlphaB)*(1-z)^(BetaB-1)/Beta(AlphaB, BetaB):
Prior := MNormDist*Dist*Dist1:

meanG1 := y*(X[1]+X[2])+X[3]+X[4]:
mu[G1] := ln(meanG1)+(1/2)*sigma[G]^2:
z[G1] := 20000;

LikG1 := exp(-(1/2)*(ln(z[G1])-mu[G1])^2/sigma[G]^2)/(sigma[G]*z[G1]*sqrt(2*Pi)):

PostNumG1 := Prior*LikG1:
DenomG1 := evalf(Int(PostNumG1, [X[1] = -infinity .. infinity, X[2] = -infinity .. infinity, X[3] = -infinity .. infinity, X[4] = -infinity .. infinity, X[5] = -infinity .. infinity, X[6] = -infinity .. infinity, X[7] = -infinity .. infinity, X[8] = -infinity .. infinity, y = 0 .. 1, z = 0 .. 1]));
evalf(DenomG1);

Thanks for any help on this!!!

Sarah

 

Please Wait...