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
TransformError::InvalidLambdaRange—lambda_min >= lambda_max, or either bound is not finiteTransformError::InsufficientData— fewer than 2 data pointsTransformError::NonPositiveData— any element ≤ 0TransformError::NonFiniteData— any element is NaN or infinite
§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);