Dear all,

I would like to use a function that computes the multivariate normal probability a in the paper "Numerical Computation of Multivariate Normal Probabilities" from Alan Genz.

This is already implemented in R, however, I want to use the Fortran function that is called from R in Octave and I can't. The closest approach was to try

mkoctfile sadmvnt.f

but the sadmvn function does not work in Octave. Can someone help me?

Simple example

Inputs en R (it works):

library("mnormt")

sadmvn(lower=c(-1.57, -0.9295, -0.0708),

upper=c(0.63,1.2705,2.1292),

mkean = c(0,0,0),

varcov = matrix(c(2.7746,-0.0757,0.0358,-0.0757,3.5240,1.3156,0.0358,1.3156,3.7585),3,3) )

Required Outputs: 0.08349516 , 1.81793e-07

R package function:

sadmvn <- function(lower, upper, mean, varcov,

maxpts=2000*d, abseps=1e-6, releps=0)

{

if(any(lower > upper)) stop("lower>upper integration limits")

if(any(lower == upper)) return(0)

d <- as.integer(if(is.matrix(varcov)) ncol(varcov) else 1)

varcov <- matrix(varcov, d, d)

sd <- sqrt(diag(varcov))

rho <- cov2cor(varcov)

lower <- as.double((lower-mean)/sd)

upper <- as.double((upper-mean)/sd)

if(d == 1) return(pnorm(upper)-pnorm(lower))

infin <- rep(2,d)

infin <- replace(infin, (upper == Inf) & (lower > -Inf), 1)

infin <- replace(infin, (upper < Inf) & (lower == -Inf), 0)

infin <- replace(infin, (upper == Inf) & (lower == -Inf), -1)

infin <- as.integer(infin)

if(any(infin == -1)) {

if(all(infin == -1)) return(1)

k <- which(infin != -1)

d <- length(k)

lower <- lower[k]

upper <- upper[k]

if(d == 1) return(pnorm(upper) - pnorm(lower))

rho <- rho[k, k]

infin <- infin[k]

if(d == 2) return(biv.nt.prob(0, lower, upper, rep(0,2), rho))

}

lower <- replace(lower, lower == -Inf, 0)

upper <- replace(upper, upper == Inf, 0)

correl <- as.double(rho[upper.tri(rho, diag=FALSE)])

maxpts <- as.integer(maxpts)

abseps <- as.double(abseps)

releps <- as.double(releps)

error <- as.double(0)

value <- as.double(0)

inform <- as.integer(0)

** result <- .Fortran("sadmvn", d, lower, upper, infin, correl, maxpts,**

** abseps, releps, error, value, inform, PACKAGE="mnormt")**

prob <- result[[10]]

attr(prob,"error") <- result[[9]]

attr(prob,"status") <- switch(1+result[[11]],

"normal completion", "accuracy non achieved", "oversize")

return(prob)

}

Fortran inputs:

d = 3

lower = -0.5658473 -0.2637628 -0.0188373

upper = 0.2270598 0.3605278 0.5665026

infin = 2 2 2

correl = -0.02420905 0.01108602 0.36149196

maxpts = 600

abseps = 1e-06

releps = 0

error = 0

value = 0

inform = 0

I appreciate very much any suggestion.

Best,

---------------

**Aurora González Vidal Ph.D. student in Data Analytics for Energy Efficiency **
T. +34 868 88 7866 Department of Information and Communication Engineering

[hidden email] Faculty of Computer Sciences

sae.saiblogs.inf.um.es University of Murcia