expand.1par.R
## Hessian, gradient expansion functions for 1-parameter base distributions
# Useful for GLM and GLM-like models
# Reference: Theorem A.1 of http://arxiv.org/pdf/1308.0657v1.pdf
### One-parameter distributions
ll_1par_expand <- function(beta,X,y,ll_1par_base) {
# obtain base distribution derivatives
ret <- ll_1par_base(X%*%beta,y)
# expand base derivatives
f <- sum(ret$f)
g <- t(X)%*%ret$g
xtw <- 0*X
for (k in 1:ncol(X)) xtw[,k] <- X[,k]*ret$h
h <- t(xtw)%*%X
return (list(f=f,g=g,h=h))
}
## log-likelihood for base GLM distributions
# Bernoulli distribution with logit link function: y ~ dBern(1/(1+exp(-u))), y=0,1
ll_bern_logit <- function(u,y) {
eu <- exp(u)
f <- -log(1+1/eu)-(1-y)*u
g <- 1/(1+eu)-(1-y)
h <- -eu/(1+eu)^2
return (list(f=f,g=g,h=h))
}
# Poisson distribution with log link function: y ~ dPois(exp(u)), y=0,1,2,...
ll_pois_log <- function(u,y) {
eu <- exp(u)
f <- y*u-eu - lfactorial(y)
g <- y-eu
h <- -eu
return (list(f=f,g=g,h=h))
}
# exponential distribution with log link function: y ~ dExp(exp(u)), y in [0,+Inf]
ll_exp_log <- function(u,y) {
eu <- exp(u)
f <- u-y*eu
g <- 1-y*eu
h <- -y*eu
return (list(f=f,g=g,h=h))
}
# goemetric distribution with logit link function: y ~ dGeom(1/(1+exp(-u)))
ll_geom_logit <- function(u,y) {
eu <- exp(u)
f <- -(y*u+(1+y)*log(1+1/eu))
g <- -y+(1+y)/(1+eu)
h <- -(1+y)*eu/(1+eu)^2
return (list(f=f,g=g,h=h))
}