Skip to main content

estimate_lambda

Function estimate_lambda 

Source
pub fn estimate_lambda(
    y: &[f64],
    lambda_min: f64,
    lambda_max: f64,
) -> Result<LambdaEstimate, TransformError>
Expand description

Estimate the optimal Box-Cox λ via maximum likelihood.

Maximises the profile log-likelihood of the normalised transform (Box & Cox 1964, §3), where g is the geometric mean of the sample:

z(λ) = g·expm1(λ·ln(y/g)) / λ        (λ ≠ 0)
       g·ln(y/g)                      (λ = 0)
ℓ(λ) = -(n/2)·ln(Var_population(z))

Normalising is what makes the likelihood computable, not merely tidier. The unnormalised (y^λ - 1)/λ puts every value next to -1 when y^λ is small – y ≈ 25, λ = -5 gives y^λ ≈ 1e-7 – so the differences between samples, which are all the variance is made of, fall off the bottom of the mantissa. Dividing through by g^(λ-1) recentres the values on zero and keeps those differences at full precision; it also folds in the Jacobian, which is why no (λ-1)·Σ ln(yᵢ) term appears above.

The search is a golden-section search over [lambda_min, lambda_max] with up to 100 iterations (terminating early when the bracket width < 1e-6).

The result says whether the maximum lies on an end of the range (LambdaEstimate::at_bound), decided by comparing the likelihood at the ends against the interior candidate rather than by where the search happened to stop. A likelihood that is monotone across the range leaves the bracket somewhere in the middle, so asking the search where it landed answers a different question than the one a caller is asking.

§Errors

§Examples

use u_numflow::transforms::estimate_lambda;

// Exponential data is well-linearised by log (λ ≈ 0).
let y: Vec<f64> = (1..=30).map(|i| (i as f64 * 0.2_f64).exp()).collect();
let est = estimate_lambda(&y, -2.0, 2.0).unwrap();
assert!(est.lambda.abs() < 0.3, "expected lambda near 0, got {}", est.lambda);
assert!(!est.at_bound);