Skip to contents

Fit a Skellam distribution with fixed mean to estimate variance

Source: R/fit_methods.R
fit_skellam.Rd

Takes in a vector of integer observations and a vector of fixed means (predictions) and estimates the Skellam distribution variance parameter. The Skellam distribution models the difference of two Poisson random variables, making it suitable for count differences that can be negative.

With fixed mean μ = λ₁ - λ₂, the variance σ² = λ₁ + λ₂ is estimated via MLE or method of moments. The Poisson rate parameters can be recovered as:

  • λ₁ = (μ + σ²) / 2

  • λ₂ = (σ² - μ) / 2

For valid parameters, we require σ² > |μ|.

Usage

fit_skellam(x, mu, method = "mle")

Arguments

x

Vector of observed integer values (can be negative).

mu

Vector of expected values (fixed means).

method

Either "mle" for maximum likelihood estimation or "mom" for method of moments. Default is "mle".

Value

The estimated variance (σ² = λ₁ + λ₂).

See also

Other estimate_observation_error: fit_normal()

Examples

# Differences can be negative
obs <- c(-2, 1, 3, 0, -1)
pred <- c(-1.5, 0.8, 2.5, 0.2, -0.5)
variance <- fit_skellam(obs, pred)
variance
#> [1] 2.510042

# To sample from Skellam with a new mean prediction:
mu_new <- 5
lambda1 <- 0.5 * (mu_new + variance)
lambda2 <- 0.5 * (variance - mu_new)
# samples <- rpois(n, lambda1) - rpois(n, lambda2)