INLAvaan (version 0.2.4)

fit_skew_normal: Fit a skew normal distribution to log-density evaluations

Description

Fit a skew normal distribution to log-density evaluations

Usage

fit_skew_normal(x, y, threshold_log_drop = -6, temp = NA)

Value

A list with fitted parameters:

• xi: location parameter

• omega: scale parameter

• alpha: shape parameter

• logC: log-normalization constant

• k: temperature parameter

• rsq: R-squared of the fit

Note that logC and k are not used when fitting from a sample.

Arguments

x

A numeric vector of points where the density is evaluated.

y

A numeric vector of log-density evaluations at points x.

threshold_log_drop

A negative numeric value indicating the log-density drop threshold below which points are ignored in the fitting. Default is -6.

temp

A numeric value for the temperature parameter k. If NA (default), it is included in the optimisation.

Details

This skew normal fitting function uses a weighted least squares approach to fit the log-density evaluations provided in y at points x. The weights are set to be the density evaluations raised to the power of the temperature parameter k. This has somewhat an interpretation of finding the skew normal fit that minimises the Kullback-Leibler divergence from the true density to it.

In R-INLA, the C code implementation from which this was translated from can be found here.

Examples

# Fit a SN curve to gamma log-density
logdens <- function(x) dgamma(x, shape = 3, rate = 1, log = TRUE)

x_grid <- seq(0.1, 8, length.out = 21)
y_log  <- sapply(x_grid, logdens)
y_log  <- y_log - max(y_log)  # normalise to have maximum at zero

res <- fit_skew_normal(x_grid, y_log, temp = 10)
unlist(res)

# Compare truth vs skew-normal approximation
x_fine <- seq(0.1, 8, length.out = 200)
y_true <- exp(logdens(x_fine))
y_sn   <- dsnorm(x_fine, xi = res$xi, omega = res$omega, alpha = res$alpha)

plot(x_fine, y_true, type = "n", xlab = "x", ylab = "Density", bty = "n")
polygon(c(x_fine, rev(x_fine)), c(y_true, rep(0, length(x_fine))),
        col = adjustcolor("#131516", 0.25), border = NA)
lines(x_fine, y_sn, col = "#00A6AA", lwd = 2)
legend("topright", legend = c("Truth", "SN Approx."),
       fill = c(adjustcolor("#131516", 0.25), NA), border = NA,
       col = c(NA, "#00A6AA"), lwd = c(NA, 2), lty = c(NA, 1),
       bty = "n")