Nothing
################################
#### Saddlepoint approximations of the Fisher-Bingham distributions
#### Tsagris Michail 02/2014
#### mtsagris@yahoo.gr
#### References: Kume Alfred and Wood Andrew T.A. (2005)
#### Saddlepoint approximations for the Bingham and Fisher-Bingham normalizing constants (Biometrika)
################################
fb.saddle <- function(gam, lam) {
## gam is the parameters of the Fisher part
## lam is the eigenvalues of the matrix of the Bingham part
lam <- sort(lam) ## sorts the eigenvalues of the Bingham part
mina <- min(lam)
if (mina <= 0) {
aaa <- abs(mina) + 1
lam <- lam + aaa ## makes all the lambdas positive and greater than zero
}
p <- length(gam) ## dimensionality of the distribution
para <- c(gam, lam) ## the parameters of the Fisher-Bingham
saddle.equat <- function(ta, para) {
## saddlepoint equation
p <- length(para)/2
gam <- para[1:p]
lam <- para[ -c(1:p) ]
sum( 0.5/(lam - ta) + 0.25 * ( gam^2/(lam - ta)^2 ) ) - 1
}
low <- lam[1] - 0.25 * p - 0.5 * sqrt(0.25 * p^2 + p * max(gam)^2) ## lower bound
up <- lam[1] - 0.25 - 0.5 * sqrt( 0.25 + min(gam)^2 ) ## not the exact upper
## bound but a bit higher
ela <- uniroot(saddle.equat, c(low, up), para = para, tol = 1e-08)
tau <- ela$root ## tau which solves the saddlepoint equation
### below are the derivatives of the cumulant generating function
kfb <- function(j, gam, lam, ta) {
if (j == 1) {
kd <- sum( 0.5/(lam - ta) + 0.25 * ( gam^2/(lam - ta)^2 ) )
} else if (j > 1) kd <- sum( 0.5 * factorial(j - 1)/(lam - ta)^j + 0.25 * factorial(j) * gam^2/(lam - ta)^(j + 1) )
kd
}
rho3 <- kfb(3, gam, lam, tau)/kfb(2, gam, lam, tau)^1.5
rho4 <- kfb(4, gam, lam, tau)/kfb(2, gam, lam, tau)^2
Ta <- rho4/8 - 5/24 * rho3^2
c1 <- 0.5 * log(2) + 0.5 * (p - 1) * log(pi) - 0.5 * log( kfb(2, gam, lam, tau) ) -
0.5 * sum( log(lam - tau) ) - tau + 0.25 * sum( gam^2/(lam - tau) )
## c1 <- sqrt(2) * pi^(0.5 * (p - 1) ) * kfb(2, gam, lam, tau)^(-0.5) *
## prod(lam - tau)^(-0.5) * exp( -tau + 0.25 * sum( gam^2/(lam - tau) ) )
c2 <- c1 + log1p(Ta)
c3 <- c1 + Ta
## the next multiplications brings the modification with the negative
## values in the lambdas back
if ( mina <= 0 ) {
c1 <- c1 + aaa
c2 <- c2 + aaa
c3 <- c3 + aaa
}
logcon <- c(c1, c2, c3)
names(logcon) <- c("first order", "second order", "third order")
logcon
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.