prior_loglik: Log-Prior Distribution Evaluation for lgspline Models

Description

Evaluates the log-prior distribution on beta coefficients conditional upon dispersion and penalaties,

Usage

prior_loglik(model_fit, sigmasq = NULL)

Value

A numeric scalar for the prior-loglikelihood (the penalty on beta coefficients actually computed)

Arguments

model_fit: An lgspline model object
sigmasq: A scalar numeric representing the dispersion parameter. By default it is NULL, and the sigmasq_tilde associated with model_fit will be used. Otherwise, custom values can be supplied.

Details

Returns the quadratic form of B^T(Lambda)B evaluated at the tuned or fixed penalties, scaled by negative one-half inverse dispersion.

Assuming fixed penalties, the prior distribution of $\beta$ is given as follows:

$$\beta | \sigma^2 \sim \mathcal{N}(\textbf{0}, \frac{1}{\sigma^2}\Lambda)$$

The log-likelihood obtained from this can be shown to be equivalent to the following, with $C$ a constant with respect to $\beta$.

$$\implies \log P(\beta|\sigma^2) = C-\frac{1}{2\sigma^2}\beta^{T}\Lambda\beta$$

This is useful for computing joint log-likelihoods and performing valid likelihood ratio tests between nested lgspline models.

Examples

Run this code

## Data
t <- sort(runif(100, -5, 5))
y <- sin(t) - 0.1*t^2 + rnorm(100)

## Model keeping penalties fixed
model_fit <- lgspline(t, y, opt = FALSE)

## Full joint log-likelihood, conditional upon known sigma^2 = 1
jntloglik <- sum(dnorm(model_fit$y,
                    model_fit$ytilde,
                    1,
                    log = TRUE)) +
          prior_loglik(model_fit, sigmasq = 1)
print(jntloglik)