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