Calculations for the Collapsed Pibble Model — loglikPibbleCollapsed • stray

Functions providing access to the Log Likelihood, Gradient, and Hessian of the collapsed pibble model. Note: These are convenience functions but are not as optimized as direct coding of the PibbleCollapsed C++ class due to a lack of Memoization. By contrast function optimPibbleCollapsed is much more optimized and massively cuts down on repeated calculations. A more efficient Rcpp module based implementation of these functions may following if the future. For model details see optimPibbleCollapsed documentation

loglikPibbleCollapsed(Y, upsilon, ThetaX, KInv, AInv, eta, sylv = FALSE)

gradPibbleCollapsed(Y, upsilon, ThetaX, KInv, AInv, eta, sylv = FALSE)

hessPibbleCollapsed(Y, upsilon, ThetaX, KInv, AInv, eta, sylv = FALSE)

Arguments

Y	D x N matrix of counts
upsilon	(must be > D)
ThetaX	D-1 x N matrix formed by Theta*X (Theta is Prior mean for regression coefficients)
KInv	Inverse of K for LTP (for Pibble defined as KInv = solve(Xi))
AInv	Inverse of A for LTP (for Pibble defined as AInv = solve(diag(N)+ t(X)%%Gamma%%X) )
eta	matrix (D-1)xN of parameter values at which to calculate quantities
sylv	(default:false) if true and if N < D-1 will use sylvester determinant identity to speed computation

Value

see below

loglikPibbleCollapsed - double
gradPibbleCollapsed - vector
hessPibbleCollapsed- matrix

Examples

D <- 10
Q <- 2
N <- 30

# Simulate Data
Sigma <- diag(sample(1:8, D-1, replace=TRUE))
Sigma[2, 3] <- Sigma[3,2] <- -1
Gamma <- diag(sqrt(rnorm(Q)^2))
Theta <- matrix(0, D-1, Q)
Phi <-  Theta + t(chol(Sigma))%*%matrix(rnorm(Q*(D-1)), nrow=D-1)%*%chol(Gamma)
X <- matrix(rnorm(N*(Q-1)), Q-1, N)
X <- rbind(1, X)
Eta <- Phi%*%X + t(chol(Sigma))%*%matrix(rnorm(N*(D-1)), nrow=D-1)
Pi <- t(driver::alrInv(t(Eta)))
Y <- matrix(0, D, N)
for (i in 1:N) Y[,i] <- rmultinom(1, sample(5000:10000), prob = Pi[,i])

# Priors
upsilon <- D+10
Xi <- Sigma*(upsilon-D)

# Precompute
KInv <- solve(Xi)
AInv <- solve(diag(N)+ t(X)%*%Gamma%*%X)
ThetaX <- Theta%*%X


loglikPibbleCollapsed(Y, upsilon, ThetaX, KInv, AInv, Eta)
#> [1] -270563.4
gradPibbleCollapsed(Y, upsilon, ThetaX, KInv, AInv, Eta)[1:5]
#> [1]   1.941015   4.862026  12.997489   6.558967 -20.942907
hessPibbleCollapsed(Y, upsilon, ThetaX, KInv, AInv, Eta)[1:5,1:5]
#>               [,1]         [,2]         [,3]          [,4]         [,5]
#> [1,] -23.447299725 -0.004375002    1.3595441    0.43590688     4.681352
#> [2,]  -0.004375002 -6.139478074    0.2997245    0.09627549     1.179519
#> [3,]   1.359544139  0.299724528 -360.2624167    6.35450008    77.884467
#> [4,]   0.435906884  0.096275494    6.3545001 -105.77379308    21.893175
#> [5,]   4.681351597  1.179519207   77.8844671   21.89317482 -1061.657850