Integration with algebraic.mle and algebraic.dist

The Three-Package Ecosystem

The likelihood.model package is part of an ecosystem of three interoperable packages:

When you fit a model with likelihood.model, the result is a fisher_mle object that inherits from algebraic.mle::mle. This means all of the algebraic.mle generics work on it directly:

library(likelihood.model)
library(algebraic.mle)

Basic Interface: params, nparams, observed_fim

Let’s fit an exponential model and explore the MLE through the algebraic.mle interface:

set.seed(42)
true_lambda <- 2.0
n <- 200
df <- data.frame(t = rexp(n, rate = true_lambda))

model <- exponential_lifetime("t")
mle_result <- fit(model)(df)

# algebraic.mle generics -- these work because fisher_mle inherits from mle
params(mle_result)
#> lambda 
#>  1.727
nparams(mle_result)
#> [1] 1

The MLE estimate is close to the true value (2.0), demonstrating good recovery with \(n = 200\) observations. The params() and nparams() generics work identically to how they work on native algebraic.mle::mle objects – because fisher_mle inherits from that class, the entire algebraic.mle interface is available without any adapter code.

The observed_fim() function returns the observed Fisher information matrix, computed as the negative Hessian of the log-likelihood at the MLE:

observed_fim(mle_result)
#>        lambda
#> lambda  67.03

For the exponential model, this is a 1x1 matrix equal to \(n/\hat\lambda^2\). It should be positive at the MLE:

cat("Observed FIM:", observed_fim(mle_result)[1,1], "\n")
#> Observed FIM: 67.03
cat("Positive:", observed_fim(mle_result)[1,1] > 0, "\n")
#> Positive: TRUE

For comparison, vcov() is the inverse of the observed FIM:

vcov(mle_result)
#>         lambda
#> lambda 0.01492
1 / observed_fim(mle_result)
#>         lambda
#> lambda 0.01492

The numerical equality confirms that vcov() is computed as solve(observed_fim()), the classical relationship \(\widehat{\text{Var}}(\hat\theta) = I(\hat\theta)^{-1}\).

Parameter Transformations via rmap

One powerful feature of algebraic.mle is the rmap() function, which transforms MLE parameters using the delta method. This propagates uncertainty through nonlinear transformations.

For an Exponential(\(\lambda\)) distribution, the mean lifetime is \(1/\lambda\). We can transform our rate MLE to get an MLE for the mean lifetime:

# Transform rate -> mean lifetime
mean_life_mle <- rmap(mle_result, function(p) {
  c(mean_lifetime = 1 / p[1])
})

params(mean_life_mle)
#> mean_lifetime 
#>        0.5789
se(mean_life_mle)
#> [1] 0.04385

Compare to the true mean lifetime:

true_mean <- 1 / true_lambda
cat("True mean lifetime:", true_mean, "\n")
#> True mean lifetime: 0.5
cat("Estimated mean lifetime:", params(mean_life_mle), "\n")
#> Estimated mean lifetime: 0.5789
cat("95% CI:", confint(mean_life_mle), "\n")
#> 95% CI: 0.493 0.6649

We can also transform to multiple derived quantities at once:

# Derive mean, variance, and median of the exponential distribution
derived_mle <- rmap(mle_result, function(p) {
  lam <- p[1]
  c(
    mean   = 1 / lam,
    var    = 1 / lam^2,
    median = log(2) / lam
  )
})

params(derived_mle)
#>   mean    var median 
#> 0.5789 0.3352 0.4013
se(derived_mle)
#> [1] 0.04137 0.04866 0.02868

Notice that the variance has the largest SE relative to its estimate – second moments are inherently harder to estimate than first moments.

Monte Carlo Expectations

The expectation() function computes Monte Carlo expectations over the MLE’s asymptotic sampling distribution.

set.seed(123)
# E[lambda^2] under the asymptotic distribution
e_lam_sq <- expectation(mle_result, function(p) p[1]^2,
                         control = list(n = 10000L))
cat("E[lambda^2]:", e_lam_sq, "\n")
#> E[lambda^2]: 2.998
cat("lambda^2 at MLE:", params(mle_result)[1]^2, "\n")
#> lambda^2 at MLE: 2.984

# Probability that lambda > 1.5 (under asymptotic distribution)
pr_lam_gt <- expectation(mle_result, function(p) as.numeric(p[1] > 1.5),
                          control = list(n = 10000L))
cat("P(lambda > 1.5):", pr_lam_gt, "\n")
#> P(lambda > 1.5): 0.9679

Mean Squared Error

Under standard MLE asymptotics, bias is zero and MSE equals the variance-covariance matrix:

mse(mle_result)
#>         lambda
#> lambda 0.01492
all.equal(mse(mle_result), vcov(mle_result))
#> [1] TRUE

The TRUE result confirms that MSE = Vcov under the assumption of zero asymptotic bias.

Bootstrap Integration

The sampler() generic creates a bootstrap sampling distribution. The resulting fisher_boot object also satisfies the mle interface:

set.seed(42)
boot_sampler <- sampler(model, df = df, par = c(lambda = 2))
boot_result <- boot_sampler(n = 200)

# Same algebraic.mle generics work
params(boot_result)
#> lambda 
#>  1.727
nparams(boot_result)
#> [1] 1
se(boot_result)
#> [1] 0.1387
bias(boot_result)
#>   lambda 
#> 0.007634

The bootstrap bias is negligible compared to the standard error, consistent with the MLE being asymptotically unbiased.

Compare bootstrap and asymptotic confidence intervals:

cat("Asymptotic 95% CI:\n")
#> Asymptotic 95% CI:
confint(mle_result)
#>         2.5% 97.5%
#> lambda 1.488 1.967

cat("\nBootstrap percentile 95% CI:\n")
#> 
#> Bootstrap percentile 95% CI:
confint(boot_result, type = "perc")
#>         2.5% 97.5%
#> lambda 1.466 2.036

Using algebraic.dist Distributions

library(algebraic.dist)

The algebraic.dist package provides distribution objects with a consistent interface. We can compare the MLE’s asymptotic sampling distribution to theoretical distribution objects.

For example, the MLE of an exponential rate parameter \(\hat\lambda\) is asymptotically normal with variance \(\lambda^2/n\). We can create this theoretical distribution:

set.seed(42)

cat("MLE:", params(mle_result), "\n")
#> MLE: 1.727
cat("SE:", se(mle_result), "\n")
#> SE: 0.1221

# Theoretical asymptotic distribution of the MLE
asymp_var <- true_lambda^2 / n
asymp_dist <- normal(mu = true_lambda, var = asymp_var)
cat("\nTheoretical asymptotic distribution:\n")
#> 
#> Theoretical asymptotic distribution:
cat("  Mean:", params(asymp_dist)[1], "\n")
#>   Mean: 2
cat("  Variance:", params(asymp_dist)[2], "\n")
#>   Variance: 0.02

# Compare: sample from the MLE's estimated distribution
mle_sampler <- sampler(mle_result)
set.seed(1)
mle_samples <- mle_sampler(5000)

# vs. sample from the theoretical distribution
dist_sampler <- sampler(asymp_dist)
set.seed(1)
dist_samples <- dist_sampler(5000)

cat("\nMLE sampler mean:", mean(mle_samples), "\n")
#> 
#> MLE sampler mean: 1.727
cat("Theoretical sampler mean:", mean(dist_samples), "\n")
#> Theoretical sampler mean: 2
cat("MLE sampler sd:", sd(mle_samples), "\n")
#> MLE sampler sd: 0.1254
cat("Theoretical sampler sd:", sd(dist_samples), "\n")
#> Theoretical sampler sd: 0.1452

Summary

The fisher_mle objects returned by likelihood.model are fully compatible with the algebraic.mle interface:

Generic What it does on fisher_mle
params() Parameter estimates (same as coef())
nparams() Number of parameters
se() Standard errors
vcov() Variance-covariance matrix
observed_fim() Observed Fisher information (\(-H\))
bias() Asymptotic bias (zero)
mse() Mean squared error (equals vcov asymptotically)
AIC() Akaike information criterion
confint() Wald confidence intervals
rmap() Delta method parameter transformations
marginal() Marginal sampling distributions
expectation() Monte Carlo expectations over sampling distribution
sampler() Asymptotic or bootstrap sampling function