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?
|
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 <[hidden email]> 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. > > ______________________________________________ > [hidden email] 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. ______________________________________________ [hidden email] 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. |
Thanks for your advice. I actually meant to ask about the "pmvt" for the distribution function. Viewing the source code "pmvt" uses the function "mvt" which uses the function "probval" which sources the fortran code:
Fortran("mvtdst", N = as.integer(n), NU = as.integer(df), LOWER = as.double(lower), UPPER = as.double(upper), INFIN = as.integer(infin), CORREL = as.double(corrF), DELTA = as.double(delta), MAXPTS = as.integer(x$maxpts), ABSEPS = as.double(x$abseps), RELEPS = as.double(x$releps), error = as.double(error), value = as.double(value), inform = as.integer(inform), PACKAGE = "mvtnorm") I wish to look at how this "mvtdst" calculates the hypergeometric function (2_F_1). Anyway that I can see that? Thanks |
On 19-03-2012, at 16:54, statfan wrote: > Thanks for your advice. I actually meant to ask about the "pmvt" for the > distribution function. Viewing the source code "pmvt" uses the function > "mvt" which uses the function "probval" which sources the fortran code: > No it doesn't "source". It call a compiled Fortran subroutine. > Fortran("mvtdst", N = as.integer(n), NU = as.integer(df), > LOWER = as.double(lower), UPPER = as.double(upper), INFIN = > as.integer(infin), > CORREL = as.double(corrF), DELTA = as.double(delta), > MAXPTS = as.integer(x$maxpts), ABSEPS = as.double(x$abseps), > RELEPS = as.double(x$releps), error = as.double(error), > value = as.double(value), inform = as.integer(inform), > PACKAGE = "mvtnorm") > > I wish to look at how this "mvtdst" calculates the hypergeometric function > (2_F_1). Anyway that I can see that? Yes. Download the source code of the package. Obtainable from CRAN. Unpack and browse. Berend ______________________________________________ [hidden email] 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. |
Free forum by Nabble | Edit this page |