Revision 5689c0b3dfd2b7ff4de58e51902745cb22cca4d5 authored by Torsten Hothorn on 08 August 1977, 00:00:00 UTC, committed by Gabor Csardi on 08 August 1977, 00:00:00 UTC
1 parent 6392497
mvt.R
pmvt <- function(lower, upper, df, corr, delta, maxpts = 25000,
abseps = 0.001, releps = 0)
{
if (df < 1)
stop("cannot compute multivariate t distribution with df < 1")
if (is.null(ncol(corr)) || length(lower) == 1)
return(list(value = pt(upper, df=df, ncp=delta) -
pt(lower, df=df, ncp=delta),
error = 0, msg="univariate: using pt"))
return(mvt(lower, upper, df, corr, delta, maxpts, abseps,releps))
}
pmvnorm <- function(lower, upper, mean, corr, maxpts = 25000,
abseps = 0.001, releps = 0)
{
if (length(mean) == 2) {
delta <- c(0,0) # bug in the Fortran Sources FIXME!
lower <- lower - mean
upper <- upper - mean
} else
delta <- mean
if (length(mean) != length(lower)) stop("wrong dimensions")
if (is.null(ncol(corr)) || length(mean) == 1)
return(list(value = pnorm(upper, mean=mean, sd=corr) -
pnorm(lower, mean=mean, sd=corr),
error = 0, msg="univariate: using pnorm"))
return(mvt(lower, upper, df=0, corr, delta, maxpts, abseps,releps))
}
mvt <- function(lower, upper, df, corr, delta, maxpts = 25000,
abseps = 0.001, releps = 0)
{
n <- ncol(corr)
if (is.null(n) || n < 2) stop("dimension less then n = 2")
if (length(lower) != n) stop("wrong dimensions")
if (length(upper) != n) stop("wrong dimensions")
if (n > 1000) stop("only dimensions 1 <= n <= 100 allowed")
infin <- rep(2, n)
infin[upper == Inf] <- 1
infin[lower == -Inf] <- 0
infin[lower == -Inf & upper == Inf] <- -1
if (n > 1) {
corrF <- matrix(as.vector(corr), ncol=n, byrow=T)
corrF <- corrF[upper.tri(corrF)]
} else corrF <- corr
lower[lower == -Inf] <- 0
upper[upper == Inf] <- 0
error <- 0; value <- 0; inform <- 0
ret <- .Fortran("mvtdst", as.integer(n), as.integer(df),
as.double(lower), as.double(upper), as.integer(infin),
as.double(corrF), as.double(delta), as.integer(maxpts),
as.double(abseps), as.double(releps),
error = as.double(error), value = as.double(value),
inform = as.integer(inform))
error <- ret$error; value <- ret$value; inform <- ret$inform
msg <- NULL
if (inform == 0) msg <- "Normal Completion"
if (inform == 1) msg <- "Completion with error > abseps"
if (inform == 3) msg <- "Covariance matrix not positive semidefinite"
if (is.null(msg)) msg <- inform
out <- list(value = value, error = error, msg = msg)
return(out)
}
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...