fortran from octave

[AURORA GONZALEZ VIDAL]

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
address@hidden                  Faculty of Computer Sciences
sae.saiblogs.inf.um.es                        University of Murcia

sadmvnt.f
Description: Text document