[R] hypergeometric function in ‘ mvtnorm’

R. Michael Weylandt michael.weylandt at gmail.com
Mon Mar 19 15:50:29 CET 2012


To view the source of (most) functions, simply type funcname without
parentheses: here, you get

> dmvt

function (x, delta, sigma, df = 1, log = TRUE, type = "shifted")
{
    if (df == 0)
        return(dmvnorm(x, mean = delta, sigma = sigma, log = log))
    if (is.vector(x)) {
        x <- matrix(x, ncol = length(x))
    }
    if (missing(delta)) {
        delta <- rep(0, length = ncol(x))
    }
    if (missing(sigma)) {
        sigma <- diag(ncol(x))
    }
    if (NCOL(x) != NCOL(sigma)) {
        stop("x and sigma have non-conforming size")
    }
    if (!isSymmetric(sigma, tol = sqrt(.Machine$double.eps),
        check.attributes = FALSE)) {
        stop("sigma must be a symmetric matrix")
    }
    if (length(delta) != NROW(sigma)) {
        stop("mean and sigma have non-conforming size")
    }
    m <- NCOL(sigma)
    distval <- mahalanobis(x, center = delta, cov = sigma)
    logdet <- sum(log(eigen(sigma, symmetric = TRUE, only.values =
TRUE)$values))
    logretval <- lgamma((m + df)/2) - (lgamma(df/2) + 0.5 * (logdet +
        m * logb(pi * df))) - 0.5 * (df + m) * logb(1 + distval/df)
    if (log)
        return(logretval)
    return(exp(logretval))
}


Most of the functions in here you can see code for the same way: the
only ones you won't be able to are eigen, lgamma, log, exp, but these
methods are pretty well-documented and you shouldn't need to find code
for them. If you do, you'll need to read the underlying C.

Michael

On Sun, Mar 18, 2012 at 11:12 PM, statfan <irene_vrbik at hotmail.com> wrote:
> Is there any way to know how the "dmvt" function computes the hypergeometric
> function needed in the calculation for the density of multivariate t
> distribution?
>
> --
> View this message in context: http://r.789695.n4.nabble.com/hypergeometric-function-in-mvtnorm-tp4483730p4483730.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.



More information about the R-help mailing list