Nothing
get_solve_kappa <- function(kappa) {
kappa_solvers <-
c("Banerjee_et_al_2005",
"Tanabe_et_al_2007",
"Sra_2012",
"Song_et_al_2012",
"uniroot",
"Newton",
"Halley",
"hybrid",
"Newton_Fourier")
if(is.numeric(kappa)) {
## CASE A: Use a common given value of kappa.
if (length(kappa) > 1)
warning("only the first element of 'kappa' is used for a common given kappa")
kappa_given <- kappa[1]
use_common_kappa <- TRUE
do_kappa <- function(norms, w, d, nu = 0)
rep.int(kappa_given, length(w))
df_kappa <- function(k) 0L
} else {
kappa <- as.list(kappa)
## Should a common value be used or not?
pos <- match("common", names(kappa), nomatch = 0L)
if(pos > 0L) {
use_common_kappa <- identical(kappa[[pos]], TRUE)
kappa <- kappa[-pos]
} else {
use_common_kappa <- FALSE
}
if(length(kappa)) {
## Solver specifications.
kname <- kappa[[1L]]
kargs <- kappa[-1L]
pos <- pmatch(tolower(kname), tolower(kappa_solvers))
if(is.na(pos))
stop("Invalid kappa solver.")
kappa <- kappa_solvers[pos]
kfun <- get(sprintf("solve_kappa_%s", kappa))
solve_kappa <- function(Rbar, d)
do.call(kfun, c(list(Rbar, d), kargs))
} else {
## Default solver.
solve_kappa <- function(Rbar, d)
solve_kappa_Newton_Fourier(Rbar, d)
}
if(use_common_kappa) {
do_kappa <- function(norms, w, d, nu = 0)
rep.int(solve_kappa(sum(norms) / (sum(w) + nu), d), length(w))
df_kappa <- function(k) 1L
} else {
do_kappa <- function(norms, w, d, nu = 0)
solve_kappa(norms / (w + nu), d)
df_kappa <- function(k) k
}
}
list(do_kappa = do_kappa,
df_kappa = df_kappa,
use_common_kappa = use_common_kappa)
}
## Try to estimate movMF via EM.
movMF <-
function(x, k, control = list(), ...)
{
## Be nice a la David.
control <- c(control, list(...))
if(missing(k)) k <- NULL
## Normalize data just in case.
x <- skmeans:::row_normalize(x)
n <- nrow(x)
d <- ncol(x)
s <- d / 2 - 1
## Control parameters.
maxiter <- control$maxiter
if(is.null(maxiter))
maxiter <- 100L
reltol <- control$reltol
if(is.null(reltol))
reltol <- sqrt(.Machine$double.eps)
verbose <- control$verbose
if(is.null(verbose))
verbose <- getOption("verbose")
E_methods <-
c("softmax",
"hardmax",
"stochmax")
E <- control$E
if(is.null(E))
E <- "softmax"
else {
pos <- pmatch(tolower(E), tolower(E_methods))
if(is.na(pos))
stop("Invalid E-step method.")
E <- E_methods[pos]
}
do_P <- switch(EXPR = E,
"softmax" =
function(G, g) exp(G - g),
"hardmax" =
function(G, g)
.posterior_from_ids(max.col(G), ncol(G)),
function(G, g)
t(apply(exp(G - g), 1,
function(prob) rmultinom(1, 1, prob))))
nu <- if (is.null(control$nu)) 0 else control$nu
solve_kappa <- get_solve_kappa(control$kappa)
do_kappa <- solve_kappa$do_kappa
df_kappa <- solve_kappa$df_kappa
use_common_kappa <- solve_kappa$use_common_kappa
if (use_common_kappa) {
if (length(nu) > 1)
warning("only the first element of 'nu' is used for common kappa")
nu <- nu[1]
}
## Allow to specify "known" ids for fitting the respective movMF
## model by maximum likelihood. This can be accomplished within our
## framework as follows:
## * Compute P as the corresponding membership matrix.
## * Perform one M step to estimate the parameters.
## * Perform one E step to compute the likelihood, but keep P.
ids <- control$ids
if(!is.null(ids)) {
if(identical(ids, TRUE))
ids <- attr(x, "z") # Be nice for the case of data
# simulated by rmovMF().
if (nrow(x) != length(ids))
stop("length of 'ids' needs to match the number of observations")
P0 <- .posterior_from_ids(ids, k)
do_P <- function(G, g) G
start <- list(P0)
maxiter <- 1L
if(!is.null(control$nruns))
warning("control argument 'nruns' ignored because ids are specified")
if(!is.null(control$start))
warning("control argument 'start' ignored because ids are specified")
} else {
## Initialization.
start <- control$start
nruns <- control$nruns
if(is.null(start)) {
if(is.null(nruns))
nruns <- 1L
start <- as.list(rep.int("p", nruns))
}
else {
if(!is.list(start))
start <- list(start)
if(!is.null(nruns))
warning("control argument 'nruns' ignored because 'start' is specified")
}
}
nruns <- length(start)
minalpha <- control$minalpha
if(is.null(minalpha)) minalpha <- 0
if(minalpha >= 1)
minalpha <- minalpha / n
converge <- control$converge
if(is.null(converge)) {
converge <- if(E == "stochmax") FALSE else TRUE
}
L_opt <- -Inf
opt_old <- opt <- NULL
run <- 1L
logLiks <- vector(length = nruns)
if(verbose && (nruns > 1L))
message(gettextf("Run: %d", run),
domain = NA)
repeat {
G <- NULL
P <- movMF_init(x, k, start[[run]])
if (!use_common_kappa) {
if (length(nu) > 1 & length(nu) != ncol(P))
warning("nu is changed to have length", ncol(P))
nu <- rep_len(nu, ncol(P))
}
L_old <- -Inf
iter <- 0L
logLiks[run] <- tryCatch({
while(iter < maxiter) {
## M step.
w <- colSums(P)
alpha <- w / n
while(any(alpha < minalpha)) {
if(verbose)
message("*** Removing one component ***")
nok <- which.min(alpha)
P <- P[, -nok, drop = FALSE]
if (!use_common_kappa) {
nu <- nu[-nok]
}
P <- do_P(P, log_row_sums(P))
alpha <- colMeans(P)
if(!is.null(G))
L_old <- sum(log_row_sums(G[, -nok, drop = FALSE]))
}
if(any(alpha == 0 & nu <= 0))
stop("Cannot handle empty components")
M <- skmeans:::g_crossprod(P, x)
norms <- row_norms(M)
M <- M / ifelse(norms > 0, norms, 1)
## If a cluster contains only identical observations,
## Rbar = 1.
kappa <- do_kappa(norms, w, d, nu)
## E step.
G <- cadd(skmeans:::g_tcrossprod(x, kappa * M),
log(alpha) - lH(kappa, s))
g <- log_row_sums(G)
L_new <- sum(g)
if(verbose && (iter %% verbose == 0))
message(gettextf("Iteration: %d *** L: %g",
iter, L_new),
domain = NA)
if(converge) {
if(abs(L_old - L_new) < reltol * (abs(L_old) + reltol)) {
L_old <- L_new
break
}
L_old <- L_new
} else if(L_new > L_old) {
L_old <- L_new
opt_old <- .movMF_object(kappa * M, alpha, L_old, P, iter)
}
P <- do_P(G, g)
iter <- iter + 1L
}
if(L_old > L_opt) {
opt <- if(converge)
.movMF_object(kappa * M, alpha, L_old, P, iter)
else opt_old
L_opt <- L_old
}
L_old
}, error = function(e) {
if(verbose) {
## <NOTE>
## Reporting problems is a bit of a mess in cases the
## above explicitly called stop(), in which case the
## condition call is
## doTryCatch(return(expr), name, parentenv, handler)
## Ideally, in these cases we would throw suitably
## classed conditions, and provide call/message methods
## for them.
call <- conditionCall(e)
msg <- conditionMessage(e)
s <- if(!is.null(call) &&
(substring(s <- deparse(call)[1L], 1L, 10L) !=
"doTryCatch"))
sprintf("Error in %s: %s", s, msg)
else
sprintf("Error: %s", msg)
message(gettextf("EM algorithm did not converge:\n%s",
s),
domain = NA)
## </NOTE>
}
NA
})
if(run >= nruns) break
run <- run + 1L
if(verbose)
message(gettextf("Run: %d", run),
domain = NA)
}
## Compute log-likelihood.
if(is.null(opt))
stop("EM algorithm did not converge for any run")
k <- length(alpha)
dimnames(opt$theta)<- list(seq_len(k), colnames(x))
dimnames(opt$P)<- list(rownames(x), seq_len(k))
ll <- L(x, opt$theta, opt$alpha)
## Add the "degrees of freedom" (number of (estimated) parameters in
## the model): with k the number of classes actually used,
## \mu: k * (d - 1) (as constrained to unit length),
## \kappa: 0, 1 or k depending on whether we use a common given,
## a common (but not given) kappa, or individual kappas.
## \alpha: k - 1 (as constrained to unit sum),
## for a total of
## k d - 1 + df_kappa(k)
attr(ll, "df") <- d * k - 1L + df_kappa(k)
attr(ll, "nobs") <- n
class(ll) <- "logLik"
opt$ll <- ll
opt$details <- list(reltol = reltol,
iter = c(iter = opt$iter, maxiter = maxiter),
logLiks = logLiks,
E = E,
kappa = control$kappa,
minalpha = minalpha,
converge = converge)
opt
}
## Generator.
.movMF_object <-
function(theta, alpha, L, P, iter)
{
o <- list(theta = theta, alpha = alpha, L = L, P = P, iter = iter)
class(o) <- "movMF"
o
}
## Methods.
print.movMF <-
function(x, ...)
{
cat("theta:\n")
print(x$theta)
cat("alpha:\n")
print(x$alpha)
cat("L:\n")
print(x$L)
invisible(x)
}
coef.movMF <-
function(object, ...)
object[c("theta", "alpha")]
logLik.movMF <-
function(object, newdata, ...)
{
if (missing(newdata))
return(object$ll)
else {
newdata <- skmeans:::row_normalize(newdata)
return(L(newdata, object$theta, object$alpha))
}
}
predict.movMF <-
function(object, newdata = NULL,
type = c("class_ids", "memberships"), ...)
{
type <- match.arg(type)
P <- if(is.null(newdata)) object$P else {
x <- skmeans:::row_normalize(newdata)
theta <- object$theta
alpha <- object$alpha
kappa <- row_norms(theta)
## Could maybe check dimensions.
## Same as for E step in movMF().
d <- nrow(theta)
s <- d / 2 - 1
G <- cadd(skmeans:::g_tcrossprod(x, theta),
log(alpha) - lH(kappa, s))
g <- log_row_sums(G)
exp(G - g)
}
if(type == "class_ids")
max.col(P)
else
P
}
## Initializer for normalized x.
## Try something similar to what we do for skmeans, but note that
## eventually we need a matrix of posterior probabilities P.
movMF_init <-
function(x, k, start)
{
if(is.character(start)) {
if(start %in% c("p", "S", "s")) {
## <FIXME>
## How should we turn initial centroids into a posterior?
## For now, do fuzzy classification with m = 2 and cosine
## dissimilarity.
M <- skmeans:::.skmeans_init_for_normalized_x(x, k, start)
D <- pmax(1 - skmeans:::g_tcrossprod(x, M), 0)
clue:::.memberships_from_cross_dissimilarities(D, 2)
## </FIXME>
}
else if(start == "i") {
## Initialize by choosing random class ids, and building a
## binary posterior from these.
ids <- sample.int(k, nrow(x), replace = TRUE)
.posterior_from_ids(ids, k)
}
else
stop("Invalid control option 'start'")
}
else if(!is.null(dim(start)))
start
else {
## A vector of class ids, hopefully.
ids <- match(start, sort(unique(start)))
.posterior_from_ids(ids, k)
}
}
## Normalizing constant for the von Mises Fisher distribution.
## (In case the reference density on the unit sphere is 1.)
C <-
function(kappa, d)
{
s <- d / 2 - 1
kappa^s / ((2 * pi)^(s + 1) * besselI(kappa, s))
}
## Utility function for scaling the columns of a matrix.
cmult <-
function(A, x)
A * rep.int(x, rep.int(nrow(A), ncol(A)))
cadd <-
function(A, x)
A + rep.int(x, rep.int(nrow(A), ncol(A)))
## Utility for computing sums avoiding *flows.
log_row_sums <-
function(x)
{
M <- x[cbind(seq_len(nrow(x)), max.col(x, "first"))]
M + log(rowSums(exp(x - M)))
}
## Utility functions for estimating kappa.
## <NOTE>
## Work only for 0 <= Rbar < 1.
## We could try to additionally implement solve_kappa(1, ...) = Inf.
## But the utilities are internal only and the other functions cannot
## gracefully handle kappa = Inf anyways ...
## </NOTE>
solve_kappa_Banerjee_et_al_2005 <-
function(Rbar, d)
{
if(any(Rbar >= 1))
stop("Cannot handle infinite concentration parameters")
Rbar * (d - Rbar ^ 2) / (1 - Rbar ^ 2)
}
solve_kappa_Tanabe_et_al_2007 <-
function(Rbar, d, c = 1, tol = 1e-6)
{
if(any(Rbar >= 1))
stop("Cannot handle infinite concentration parameters")
old <- Rbar * (d - c) / (1 - Rbar ^ 2)
repeat {
kappa <- ifelse(old == 0, 0, old * (Rbar / A(old, d)))
if(max(abs(kappa - old)) < tol) break
old <- kappa
}
kappa
}
solve_kappa_Sra_2012 <-
function(Rbar, d)
{
## Initialize using the Banerjee et al approximation.
kappa <- solve_kappa_Banerjee_et_al_2005(Rbar, d)
## Perform two Newton steps.
A <- A(kappa, d)
kappa <- kappa - (A - Rbar) / Aprime(kappa, d, A = A)
A <- A(kappa, d)
kappa <- kappa - (A - Rbar) / Aprime(kappa, d, A = A)
kappa
}
solve_kappa_Halley <-
function(Rbar, d, tol = 1e-6, maxiter = 100L)
{
if(any(Rbar >= 1))
stop("Cannot handle infinite concentration parameters")
n <- max(length(Rbar), length(d))
d <- rep_len(d, n)
Rbar <- rep_len(Rbar, n)
sapply(seq_along(Rbar),
function(i) {
r <- Rbar[i]
r2 <- r^2
r3 <- r2 * r
r4 <- r3 * r
D <- d[i]
Dm1 <- D - 1
nu <- D / 2 - 1
kappa_0 <- Inf
kappa <- Rinv_lower_Amos_bound(r, nu)
A <- A(kappa, D)
Adiff <- A - r
for (n in seq_len(maxiter)) {
Adiff2 <- Adiff^2
Adiff3 <- Adiff^3
Dkappa <- Dm1 / kappa
numerator <-
2 * (Adiff * (1 - r * Dkappa - r2) - Adiff2 * (2 * r + Dkappa) - Adiff3)
denominator <-
2 * (r4 + r3 * 2 * Dkappa + r2 * (-2 + Dkappa^2) - r * 2 * Dkappa + 1) +
Adiff * (6 * r3 + 9 * r2 * Dkappa + r * (-6 + Dkappa * (3 * D - 4) / kappa) - 3 * Dkappa) +
Adiff2 * (6 * r2 + r * 6 * Dkappa + Dkappa * (D - 2)/kappa - 2) +
Adiff3 * (2 * r + Dkappa)
kappa <- kappa - numerator / denominator
A <- A(kappa, D)
Adiff <- A - r
if ((abs(Adiff) < tol) ||
(abs(kappa_0 - kappa) < tol * (kappa_0 + tol))) {
break
}
kappa_0 <- kappa
}
kappa
})
}
solve_kappa_Song_et_al_2012 <-
function(Rbar, d)
{
## Initialize using the Banerjee et al approximation.
kappa <- solve_kappa_Banerjee_et_al_2005(Rbar, d)
## Perform two Halley steps.
A <- A(kappa, d)
Adiff <- A - Rbar
Aprime <- Aprime(kappa, d, A = A)
kappa <- kappa - 2 * Adiff * Aprime /
(2 * Aprime^2 - Adiff * Adoubleprime(kappa, d, A = A))
A <- A(kappa, d)
Adiff <- A - Rbar
Aprime <- Aprime(kappa, d, A = A)
kappa <- kappa - 2 * Adiff * Aprime /
(2 * Aprime^2 - Adiff * Adoubleprime(kappa, d, A = A))
kappa
}
solve_kappa_uniroot <-
function(Rbar, d, tol = 1e-6)
{
if(any(Rbar >= 1))
stop("Cannot handle infinite concentration parameters")
n <- max(length(Rbar), length(d))
d <- rep_len(d, n)
Rbar <- rep_len(Rbar, n)
sapply(seq_along(Rbar),
function(i) {
r <- Rbar[i]
D <- d[i]
nu <- D / 2 - 1
interval <- c(Rinv_lower_Amos_bound(r, nu),
Rinv_upper_Amos_bound(r, nu))
if (abs(diff(interval)) < tol)
mean(interval)
else
uniroot(function(kappa) A(kappa, D) - r,
interval = interval,
tol = tol)$root})
}
solve_kappa_Newton <-
function(Rbar, d, tol = 1e-6, maxiter = 100)
{
if(any(Rbar >= 1))
stop("Cannot handle infinite concentration parameters")
n <- max(length(Rbar), length(d))
d <- rep_len(d, n)
Rbar <- rep_len(Rbar, n)
sapply(seq_along(Rbar),
function(i) {
r <- Rbar[i]
D <- d[i]
nu <- D / 2 - 1
kappa_0 <- Inf
kappa <- Rinv_lower_Amos_bound(r, nu)
for (n in seq_len(maxiter)) {
A <- A(kappa, D)
Adiff <- A - r
kappa <- kappa - Adiff / Aprime(kappa, D, A = A)
A <- A(kappa, D)
if ((abs(Adiff) < tol) ||
(abs(kappa_0 - kappa) < tol * (kappa_0 + tol)))
break
kappa_0 <- kappa
}
kappa
})
}
solve_kappa_Newton_Fourier <-
function(Rbar, d, tol = 1e-6, maxiter = 100L)
{
if(any(Rbar >= 1))
stop("Cannot handle infinite concentration parameters")
n <- max(length(Rbar), length(d))
d <- rep_len(d, n)
Rbar <- rep_len(Rbar, n)
sapply(seq_along(Rbar),
function(i) {
r <- Rbar[i]
D <- d[i]
nu <- D / 2 - 1
lower <- Rinv_lower_Amos_bound(r, nu)
upper <- Rinv_upper_Amos_bound(r, nu)
iter <- 1L
while(iter <= maxiter) {
A <- A(lower, D)
Aprime <- Aprime(lower, D, A = A)
lower <- lower - (A - r) / Aprime
A <- A(upper, D)
upper <- upper - (A - r) / Aprime
if((upper - lower) < tol * (lower + upper)) {
if ((upper - lower) < - tol * (lower + upper))
stop("no convergence")
break
}
iter <- iter + 1L
}
## <FIXME>
## What should we really return?
(lower + upper) / 2
## </FIXME>
})
}
## Combination of Newton with bisection
## See Press et al., Numerical Recipes
## Modified to
## allow different methods for determining steps (Newton, Halley)
step_Newton <-
function(x, d, Rbar, A = NULL)
{
if (is.null(A))
A <- A(x, d)
fx <- A - Rbar
dfx <- Aprime(x, d, A = A)
fx / dfx
}
step_Halley <-
function(x, d, Rbar, A = NULL)
{
if (is.null(A))
A <- A(x, d)
Adiff <- A - Rbar
Adiff2 <- Adiff^2
Adiff3 <- Adiff2 * Adiff
dx <- (d - 1) / x
Rbar2 <- Rbar^2
Rbar3 <- Rbar2 * Rbar
Rbar4 <- Rbar3 * Rbar
numerator <- 2 * (Adiff * (1 - Rbar * dx - Rbar2) - Adiff2 * (2 * Rbar + dx) - Adiff3)
denominator <-
2 * (Rbar4 + Rbar3 * 2 * dx + Rbar2 * (-2 + dx^2) - Rbar * 2 * dx + 1) +
Adiff * (6 * Rbar3 + 9 * Rbar2 * dx + Rbar * (-6 + dx * (3 * d - 4) / x) - 3 * dx) +
Adiff2 * (6 * Rbar2 + Rbar * 6 * dx + dx * (d - 2)/x - 2) +
Adiff3 * (2 * Rbar + dx)
numerator / denominator
}
solve_kappa_hybrid <-
function(Rbar, d, tol = 1e-6, maxiter = 100, step = step_Halley)
{
if(any(Rbar >= 1))
stop("Cannot handle infinite concentration parameters")
n <- max(length(Rbar), length(d))
d <- rep_len(d, n)
Rbar <- rep_len(Rbar, n)
sapply(seq_along(Rbar),
function(i) {
r <- Rbar[i]
D <- d[i]
nu <- D / 2 - 1
rtsold <- xl <- rts <- Rinv_lower_Amos_bound(r, nu)
xh <- Rinv_upper_Amos_bound(r, nu)
dxold <- abs(xh - xl)
Arts <- A(rts, D)
for (j in seq_len(maxiter)) {
dx <- step(rts, D, r, Arts)
rts <- rts - dx
if ((rts > xh) || (rts < xl) || (abs(2 * dx) > abs(dxold))) {
dx <- 0.5 * (xh - xl)
rts <- xl + dx
if(any(xl == rts))
break
}
else {
if (rtsold == rts)
break
}
if (abs(dx) < tol)
break
dxold <- dx
rtsold <- rts
Arts <- A(rts, D)
if ((Arts - r) < 0.0)
xl <- rts
else
xh <- rts
}
rts
})
}
## Utility functions for computing A, H and its logarithm (see the
## implementation notes).
##
## With R_\nu = I_{\nu+1} / I_\nu the commonly considered modified
## Bessel function ratio and
##
## log(H_\nu)(\kappa) = \int_0^\kappa R_\nu(t) dt
##
## (so that (log(H_\nu))' = R_\nu and H_\nu(0) = 1), we have
##
## A_d = I_{d/2} / I_{d/2-1} = R_{d/2-1}
##
## and
##
## log(f(x|theta)) = theta' x - log(H_{d/2-1})(\|theta\|).
##
## For R_\nu we have Amos-type bounds
##
## G_{\alpha,\beta}(t) = t / (alpha + sqrt(t^2 + beta^2))
##
## where G_{\nu+1/2,\nu+3/2} is the uniformly best lower Amos-type
## bound, whereas there is no uniformly optimal upper bound, with
## G_{\nu,\nu+2} and G_{\nu+1/2,beta_SS(\nu)} being second-order exact
## and optimal at 0 and infinity, respectively, where
##
## beta_SS(\nu) = sqrt((\nu+1/2)(\nu+3/2)).
##
## Inverses and anti-derivatives of G are given by
##
## G_{\alpha,\beta}^{-1}(\rho)
## = \frac{\rho}{1 - \rho^2}
## (\alpha + \sqrt{\alpha^2 \rho^2 + \beta^2 (1 - \rho^2)})
##
## and
##
## S_{\alpha,\beta}(\kappa)
## = \sqrt{\kappa^2 + \beta^2} - \beta
## - \alpha \log(\alpha + \sqrt{\kappa^2 + \beta^2})
## + \alpha \log(\alpha + \beta).
##
## We use
##
## \max(G_{\nu,\nu+2}^{-1}, G_{\nu+1/2,\beta_{SS}(\nu)}^{-1})
##
## as lower bound and approximation for R_\nu^{-1},
##
## G_{\nu+1/2,\nu+3/2}^{-1}
##
## as upper bound for R_{\nu}^{-1}, and the antiderivative of
## \min(G_{\nu,\nu+2},G_{\nu+1/2,beta_SS(\nu)}) which is given by
##
## S_{\nu+1/2,\beta_{SS}(\nu)}(kappa)
## - S_{\nu+1/2,\beta_{SS}(\nu)}(\min(kappa, \kappa_\nu))
## + S_{\nu,\nu+2}(\min(kappa, \kappa_\nu))
##
## where
##
## \kappa_\nu = \sqrt{(3\nu + 11/2) (\nu + 3/2)}
##
## is the positive solution of
##
## G_{\nu,\nu+2}(\kappa) = G_{\nu+1/2,beta_SS(\nu)}(\kappa)
##
## as approximation for \log(H_\nu).
G <-
function(kappa, alpha, beta)
{
n <- max(length(kappa), length(alpha), length(beta))
kappa <- rep_len(kappa, n)
alpha <- rep_len(alpha, n)
beta <- rep_len(beta, n)
y <- numeric(n)
ind <- ((alpha == 0) & (beta == 0))
if(any(ind)) y[ind] <- 1
ind <- !ind
if(any(ind)) {
denom <- alpha[ind] + sqrt(kappa[ind]^2 + beta[ind]^2)
y[ind] <- ifelse(denom == 0, NaN, kappa[ind] / denom)
}
y
}
beta_SS <-
function(nu)
sqrt((nu + 1/2) * (nu + 3/2))
Ginv <-
function(rho, alpha, beta)
{
## Perhaps recycle arguments eventually ...
sigma <- rho^2
rho * (alpha + sqrt(alpha^2 * sigma + beta^2 * (1 - sigma))) /
(1 - sigma)
}
Rinv_lower_Amos_bound <-
function(rho, nu)
{
## Perhaps recycle arguments eventually ...
pmax(Ginv(rho, nu, nu + 2),
Ginv(rho, nu + 1/2, beta_SS(nu)))
}
Rinv_upper_Amos_bound <-
function(rho, nu)
{
## Perhaps recycle arguments eventually ...
Ginv(rho, nu + 1/2, nu + 3/2)
}
S <-
function(kappa, alpha, beta)
{
n <- max(length(kappa), length(alpha), length(beta))
kappa <- rep_len(abs(kappa), n)
alpha <- rep_len(alpha, n)
beta <- rep_len(abs(beta), n)
ab <- alpha + beta
s <- double(n)
ind <- (ab < 0) & (kappa^2 > alpha^2 - beta^2)
if (any(ind))
s[ind] <- NaN
ind <- !ind
if (any(ind)) {
u <- sqrt(kappa[ind]^2 + beta[ind]^2)
s[ind] <- u - beta[ind] - ifelse(alpha[ind] == 0, 0, alpha[ind] * log((alpha[ind] + u) / ab[ind]))
}
s
}
## Utility functions for computing H and log(H).
H <-
function(kappa, nu, v0 = 1)
{
## Compute v0 H_\nu(\kappa) by direct Taylor series summation.
## Add range checking eventually.
n <- max(length(kappa), length(nu), length(v0))
kappa <- rep_len(kappa, n)
nu <- rep_len(nu, n)
v0 <- rep_len(v0, n)
.C(C_my0F1,
as.integer(n),
as.double(kappa ^ 2 / 4),
as.double(nu + 1),
as.double(v0),
y = double(n))$y
}
lH_asymptotic <-
function(kappa, nu)
{
## Compute a suitable asymptotic approximation to
## \log(H_\nu(\kappa)).
## Add range checking eventually.
n <- max(length(kappa), length(nu))
kappa <- rep_len(kappa, n)
nu <- rep_len(nu, n)
y <- double(n)
ipk <- (kappa > 0)
ipn <- (nu > 0)
ind <- ipk & !ipn
if(any(ind)) {
## For \log(H_0) = \log(I_0), use the asymptotic approximation
## I_0(\kappa) \approx e^\kappa / \sqrt{2 \pi \kappa}
## (e.g., http://dlmf.nist.gov/10.40).
y[ind] <- kappa[ind] - log(2 * pi * kappa[ind]) / 2
}
ind <- ipk & ipn
if(any(ind)) {
## For \nu > 0, use the Amos-type based approximation discussed
## above.
kappa <- kappa[ind]
nu <- nu[ind]
beta <- beta_SS(nu)
kappa_L <- pmin(kappa, sqrt((3 * nu + 11 / 2) * (nu + 3 / 2)))
y[ind] <- S(kappa, nu + 1/2, beta) +
(S(kappa_L, nu, nu + 2) - S(kappa_L, nu + 1/2, beta))
}
y
}
lH <-
function(kappa, nu)
{
## Compute \log(H_\nu(\kappa)) (or an approximation to it).
## See the implementation notes for details.
n <- max(length(kappa), length(nu))
kappa <- rep_len(kappa, n)
nu <- rep_len(nu, n)
y <- lH_asymptotic(kappa, nu)
## If the value from the asymptotic approximation is small enough,
## we can use direct Taylor series summation.
ind <- y <= 699.5
if(any(ind))
y[ind] <- log(H(kappa[ind], nu[ind]))
## For intermediate values of the asymptotic approximation, we use
## rescaling and direct Taylor series summation.
ind <- !ind & (y <= 1399)
if(any(ind)) {
v <- y[ind] / 2
y[ind] <- v + log(H(kappa[ind], nu[ind], exp(-v)))
}
## (Otherwise, we use the asymptotic approximation.)
y
}
## For testing this: delta should be close to 0.
delta <-
function(kappa, nu)
{
(besselI(kappa, nu) -
H(kappa, nu) * (kappa / 2)^nu / gamma(nu + 1))
}
## A, A prime and A double prime.
A <-
function(kappa, d, method = c("PCF", "GCF", "RH"), tol = 1e-6)
{
method <- match.arg(method)
n <- max(length(kappa), length(d))
kappa <- rep_len(kappa, n)
d <- rep_len(d, n)
A <- vector("numeric", length = n)
index <- kappa >= tol
if (sum(index)) {
method <- match.arg(method)
if(method == "PCF") {
## Use the Perron continued fraction for R_{\nu-1}.
## Implementation based on Eqn 3.3' in Gautschi and Slavik
## (1978).
## Note that A_d = R_{d/2-1}.
A[index] <- .C(C_mycfP,
as.integer(sum(index)),
as.double(kappa[index]),
as.double(d[index] / 2),
y = double(sum(index)))$y
}
else if(method == "GCF") {
## Use the Gauss continued fraction for R_{\nu-1}.
## Implementation based on Eqn 3.2' in Gautschi and Slavik
## (1978).
## Note that A_d = R_{d/2-1}.
A[index] <- .C(C_mycfG,
as.integer(sum(index)),
as.double(kappa[index]),
as.double(d[index] / 2),
y = double(sum(index)))$y
}
else {
## Compute A_d via a ratio of H functions:
## A_d(\kappa)
## = (kappa/d) * H_{d/2}(\kappa) / H_{d/2-1}(\kappa).
s <- d[index] / 2 - 1
kappai <- kappa[index]
y <- kappai / d[index]
a <- lH_asymptotic(kappai, s + 1)
ind <- (a <= 699.5)
if(any(ind))
y[ind] <- y[ind] *
(H(kappai[ind], s + 1) /
H(kappai[ind], s))
ind <- !ind & (a <= 1399)
if(any(ind)) {
v <- exp(- a[ind] / 2)
y[ind] <- y[ind] *
(H(kappai[ind], s + 1, v) /
H(kappai[ind], s, v))
}
ind <- (a > 1399)
if(any(ind))
y[ind] <- y[ind] *
exp(a[ind] - lH_asymptotic(kappai[ind], s))
if(any(y >= 1))
stop("RH evaluation gave infeasible values which are not in the range [0, 1)")
A[index] <- y
}
}
if (sum(!index)) {
di <- d[!index]
kappai <- kappa[!index]
A[!index] <- kappai / di - kappai^3 / (di^2 * (di + 2)) + 2 * kappai^5 / (di^3 * (di + 2) * (di + 4))
}
A
}
Aprime <-
function(kappa, d, method = c("PCF", "GCF", "RH"), tol = 1e-6, A = NULL)
{
n <- max(length(kappa), length(d), length(A))
kappa <- rep_len(kappa, n)
d <- rep_len(d, n)
if(!is.null(A))
A <- rep_len(A, n)
aprime <- vector("numeric", length = n)
index <- kappa >= tol
if (sum(index)) {
if(is.null(A))
A <- A(kappa[index], d[index], method, tol)
else
A <- A[index]
aprime[index] <- 1 - A ^ 2 - A * (d[index] - 1) / kappa[index]
}
if (sum(!index)) {
di <- d[!index]
aprime[!index] <- 1 / di - 3 / (di^2 * (di + 2)) * kappa[!index]^2 + 10 / (di^3 * (di + 2) * (di + 4)) * kappa[!index]^4
}
aprime
}
Adoubleprime <-
function(kappa, d, method = c("PCF", "GCF", "RH"), tol = 1e-6, A = NULL)
{
n <- max(length(kappa), length(d), length(A))
kappa <- rep_len(kappa, n)
d <- rep_len(d, n)
if(!is.null(A))
A <- rep_len(A, n)
adoubleprime <- vector("numeric", length = n)
index <- kappa >= tol
if (sum(index)) {
di <- d[index]
if (is.null(A))
A <- A(kappa[index], di, method, tol)
else
A <- A[index]
kappa2 <- kappa[index]^2
adoubleprime[index] <- 2 * A^3 + 3 * (di - 1) / kappa[index] * A^2 + (di^2 - di - 2 * kappa2) / kappa2 * A - (di - 1) / kappa[index]
}
if(sum(!index)) {
di <- d[!index]
adoubleprime[!index] <- - 6 / (di^2 * (di + 2)) * kappa[!index] + 40 / (di^3 * (di + 2) * (di + 4)) * kappa[!index]^3
}
adoubleprime
}
## R, R prime and R double prime.
R <-
function(t, nu, ...)
A(t, 2 * (nu + 1), ...)
Rprime <-
function(t, nu, ..., R = NULL)
Aprime(t, 2 * (nu + 1), ..., A = R)
Rdoubleprime <-
function(t, nu, ..., R = NULL)
Adoubleprime(t, 2 * (nu + 1), ..., A = R)
Rtripleprime <-
function(t, nu, ..., tol = 1e-6)
{
n <- max(length(t), length(nu))
t <- rep_len(t, n)
nu <- rep_len(nu, n)
rnu_tripleprime <- vector("numeric", length = n)
ind <- t >= tol
if(sum(ind)) {
t_ind <- t[ind]
nu_ind <- nu[ind]
rnu <- R(t_ind, nu_ind, ..., tol = tol)
## No given R for now.
rnu_prime <- Rprime(t_ind, nu_ind, tol = tol, R = rnu)
rnu_doubleprime <- Rdoubleprime(t_ind, nu_ind, tol = tol, R = rnu)
delta <- 2 * nu_ind + 1
rnu_tripleprime[ind] <-
(- 2 * delta * rnu / t_ind^3
+ 2 * delta * rnu_prime / t_ind^2
- delta * rnu_doubleprime / t_ind
- 2 * (rnu_prime^2 + rnu * rnu_doubleprime))
}
ind <- !ind
if(sum(ind)) {
t_ind <- t[ind]
nu_ind <- nu[ind]
rnu_tripleprime[ind] <-
(- 0.75 / ((nu_ind + 1)^2 * (nu_ind + 2))
+ 3.75 * t_ind^2 /
( (nu_ind + 1)^3 * (nu_ind + 2) * (nu_ind + 3) ))
}
rnu_tripleprime
}
## Log-likelihood
L <-
function(x, theta, alpha)
{
sum(ldmovMF(x, theta, alpha))
}
## Utility for computing a binary posterior from the class ids.
.posterior_from_ids <-
function(ids, k = NULL)
{
if(is.null(k))
k <- max(ids)
else if(max(ids) < k) {
k <- max(ids)
warning(gettextf("due to initialization number of components reduced to %d",
k),
domain = NA)
}
else if(max(ids) > k)
stop("number of components k smaller than those provided for initialization")
n <- length(ids)
P <- matrix(0, n, k)
P[cbind(seq_len(n), ids)] <- 1
P
}
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.