fbase1.poisson.log {RegressionFactory} | R Documentation |
Single-Parameter Base Log-likelihood Function for Poisson GLM
Description
Vectorized, single-parameter base log-likelihood functions for Poisson GLM using log link function. The base function(s) can be supplied to the expander function regfac.expand.1par
in order to obtain the full, high-dimensional log-likleihood and its derivatives.
Usage
fbase1.poisson.log(u, y, fgh=2)
Arguments
u |
Varying parameter of the base log-likelihood function. This parameter is intended to be projected onto a high-dimensional space using the familiar regression transformation of |
y |
Fixed slot of the base distribution, corresponding to the response variable in the regression model. For |
fgh |
Integer with possible values 0,1,2. If |
Value
If fgh==0
, the function returns y*u-exp(u)-lfactorial(y)
for log
. If fgh==1
, a list is returned with elements f
and g
, where the latter is a vector of length length(u)
, with each element being the first derivative of the above expressions. If fgh==2
, the list will include an element named h
, consisting of the second derivatives of f
with respect to u
.
Note
In all base log-likelihood functions, we have dropped any additive terms that are independent of the distribution parameter, e.g. constant terms or those terms that are dependent on the response variable only. This is done for computational efficiency. Therefore, these functions cannot be used to obtain the absolute values of log-likelihood functions but only in the context of optimization and/or sampling. Users can write thin wrappers around these functions to add the constant terms to the function value. (Derivatives do not need correction. For Poisson family, the lfactorial(y)
term is dropped.)
Author(s)
Alireza S. Mahani, Mansour T.A. Sharabiani
See Also
Examples
## Not run:
library(sns)
library(MfUSampler)
# using the expander framework and base distributions to define
# log-likelihood function for Poisson regression
loglike.poisson <- function(beta, X, y, fgh) {
regfac.expand.1par(beta, X, y, fbase1.poisson.log, fgh)
}
# generate data for Poisson regression
N <- 1000
K <- 5
X <- matrix(runif(N*K, min=-0.5, max=+0.5), ncol=K)
beta <- runif(K, min=-0.5, max=+0.5)
y <- rpois(N, lambda = exp(X%*%beta))
# obtaining glm coefficients for comparison
beta.glm <- glm(y~X-1, family="poisson")$coefficients
# mcmc sampling of log-likelihood
nsmp <- 100
# Slice Sampler (no derivatives needed)
beta.smp <- array(NA, dim=c(nsmp,K))
beta.tmp <- rep(0,K)
for (n in 1:nsmp) {
beta.tmp <- MfU.Sample(beta.tmp
, f=loglike.poisson, X=X, y=y, fgh=0)
beta.smp[n,] <- beta.tmp
}
beta.slice <- colMeans(beta.smp[(nsmp/2+1):nsmp,])
# Adaptive Rejection Sampler
# (only first derivative needed)
beta.smp <- array(NA, dim=c(nsmp,K))
beta.tmp <- rep(0,K)
for (n in 1:nsmp) {
beta.tmp <- MfU.Sample(beta.tmp, uni.sampler="ars"
, f=function(beta, X, y, grad) {
if (grad)
loglike.poisson(beta, X, y, fgh=1)$g
else
loglike.poisson(beta, X, y, fgh=0)
}
, X=X, y=y)
beta.smp[n,] <- beta.tmp
}
beta.ars <- colMeans(beta.smp[(nsmp/2+1):nsmp,])
# SNS (Stochastic Newton Sampler)
# (both first and second derivative needed)
beta.smp <- array(NA, dim=c(nsmp,K))
beta.tmp <- rep(0,K)
for (n in 1:nsmp) {
beta.tmp <- sns(beta.tmp, fghEval=loglike.poisson, X=X, y=y, fgh=2, rnd = n>nsmp/4)
beta.smp[n,] <- beta.tmp
}
beta.sns <- colMeans(beta.smp[(nsmp/2+1):nsmp,])
# compare results
cbind(beta.glm, beta.slice, beta.ars, beta.sns)
## End(Not run)