R/diri.nr.R In Compositional: Compositional Data Analysis

Documented in diri.nr

################################
#### Dirichlet distribution parameters
#### via Newton-Raphson
#### Tsagris Michail 8/2015

#### References: Estimating a Dirichlet distribution (2012)
#### Thomas P. Minka
#### http://research.microsoft.com/en-us/um/people/minka/papers/dirichlet/minka-dirichlet.pdf
################################
diri.nr <- function(x, type = 1, tol = 1e-07) {

if (type == 1) {
runtime <- proc.time()
n <- dim(x)[1]  ## the sample size
p <- dim(x)[2]  ## dimensionality
m <- Rfast::colmeans(x)
zx <- t( log(x) )
down <-  - sum( m * ( Rfast::rowmeans( zx ) - log(m) ) )
sa <- 0.5 * (p - 1) / down  ## initial value for precision
a1 <- sa * m  ## initial values
gm <- Rfast::rowsums(zx)
z <- n * digamma( sa )
g <- z - n * digamma(a1) + gm
qk <-  - n * trigamma(a1)
b <- sum(g / qk) / ( 1/z - sum(1 / qk) )
a2 <- a1 - (g - b)/qk
i <- 2
while ( sum( abs( a2 - a1 ) ) > tol ) {
i <- i + 1
a1 <- a2
z <- n * digamma( sum(a1) )
g <- z - n * digamma(a1) + gm
qk <-  - n * trigamma(a1)
b <- sum(g / qk) / ( 1/z - sum(1 / qk) )
a2 <- a1 - (g - b) / qk
}
loglik <- n * lgamma( sum(a2) ) - n * sum( lgamma(a2) ) + sum( zx * (a2 - 1) )
runtime <- proc.time() - runtime

if ( is.null(colnames(x)) ) {
names(a2) <- paste("X", 1:p, sep = "")
} else  names(a2) <- colnames(x)
res <- list(iter = i, loglik = loglik, param = a2, runtime = runtime)

} else  res <- Rfast::diri.nr2(x, tol = tol)

res
}

Try the Compositional package in your browser

Any scripts or data that you put into this service are public.

Compositional documentation built on March 18, 2018, 1:57 p.m.