Distributional (Location-Scale): Thin Delegation

Navigation:

Implementation only, for the two-stage schedule and quantile formulas see the math theory.

There is no DistributionalTAM class: a dict formula makes StaticTAM a thin frontend that delegates the whole schedule to statistics/estimation/_distributional.py. Sub-models are built with type(model)(...), so the module never imports StaticTAM (no circular import).

The schedule (location, then scale on the residuals)

def fit(model, data: pd.DataFrame, select: str = "fixed",
        cv_alpha_grid: Sequence[float] = (-6.0, -4.0, -2.0, 0.0), cv_fraction: float = 0.3):
    """Fit the location, then the scale on its squared residuals, then the tail family.

    select chooses the location smoothing: 'fixed' uses location_alpha_p; 'gcv' runs auto_fit
    (exact for the location); 'cv' picks it from cv_alpha_grid by held-out error. The scale
    smoothing stays fixed (Gaussian GCV is invalid for the Gamma GLM).
    """
    working = data.copy()
    working["__mu__"] = _to_model_scale(model, working[model.target_col_].to_numpy())

    location_alpha = model._location_alpha_p_ if model._location_alpha_p_ is not None else model.default_alpha_p_
    if select == "gcv":
        model._location_submodel_ = _build_location_submodel(model, location_alpha)
        try:
            model._location_submodel_.auto_fit(working)
        except Exception:
            model._location_submodel_.fit(working)
    elif select == "cv":
        location_alpha = _select_location_alpha_by_cv(model, working, cv_alpha_grid, cv_fraction)
        model._location_submodel_ = _build_location_submodel(model, location_alpha).fit(working)
    else:
        model._location_submodel_ = _build_location_submodel(model, location_alpha).fit(working)

    mu_hat = _estimated(model._location_submodel_.predict(working), "__mu__")
    residual = working["__mu__"].to_numpy() - mu_hat
    squared_residual = np.clip(residual ** 2, _TINY, None)
    model._global_scale_variance_ = float(np.mean(squared_residual))

    if model._scale_loss_ == "gamma":
        working["__sigma__"] = squared_residual
    else:
        working["__sigma__"] = np.log(squared_residual + _TINY)

    model._scale_submodel_ = _build_scale_submodel(model).fit(working)

    _, sigma_hat = mu_sigma(model, working)
    fit_tail_family(model, residual / sigma_hat)
    return model

The tail law and the location/scale predictions

def select_tail_family(
    standardized_residual: np.ndarray,
    tail_family: str,
    kurtosis_threshold: float,
) -> Tuple[str, Optional[float]]:
    """Choose the standardized-residual law by excess kurtosis.

    'normal' / 'student_t' force the family; 'auto' picks Student-t only when the excess kurtosis
    exceeds kurtosis_threshold, fitting nu from it. Returns (family, nu); nu is None for the Normal.
    """
    centred = standardized_residual - standardized_residual.mean()
    variance = np.mean(centred ** 2)
    excess_kurtosis = float(np.mean(centred ** 4) / max(variance ** 2, _TINY) - 3.0)
    if tail_family == "normal":
        return "normal", None
    if tail_family == "student_t":
        return "student_t", float(np.clip(4.0 + 6.0 / max(excess_kurtosis, 1e-3), 3.0, 60.0))
    # "auto"
    if excess_kurtosis <= kurtosis_threshold:
        return "normal", None
    return "student_t", float(np.clip(4.0 + 6.0 / excess_kurtosis, 3.0, 60.0))

def mu_sigma(model, data: pd.DataFrame) -> Tuple[np.ndarray, np.ndarray]:
    """Return the location mu_hat (model scale) and scale sigma_hat for each row of data."""
    mu_hat = _estimated(model._location_submodel_.predict(data), "__mu__")
    scale_prediction = _estimated(model._scale_submodel_.predict(data), "__sigma__")
    if model._scale_loss_ == "gamma":
        conditional_variance = np.clip(scale_prediction, _TINY, None)
    else:
        conditional_variance = np.exp(scale_prediction - _LOG_CHI2_1_MEAN)
    if model.scale_shrinkage > 0.0 and model._global_scale_variance_ is not None:
        conditional_variance = (
            model.scale_shrinkage * model._global_scale_variance_
            + (1.0 - model.scale_shrinkage) * conditional_variance
        )
    return mu_hat, np.clip(np.sqrt(conditional_variance), _TINY, None)

For a log target the quantile and the median are exp of a model-scale value. A value beyond the float64 range (a location or scale prediction far outside the training range) is returned as inf with an explicit UserWarning naming the quantity and the first rows, never numpy’s silent overflow warning:

def _to_response_scale(model_scale: np.ndarray, quantity: str) -> np.ndarray:
    """exp() of a model-scale value for a log target; a value beyond float64 is inf, with an explicit warning (not numpy's)."""
    with np.errstate(over="ignore"):
        response = np.exp(model_scale)
    overflowed = np.flatnonzero(np.isposinf(response) & (model_scale > 0))
    if overflowed.size:
        warnings.warn(
            f"exp overflow in {quantity}: {overflowed.size} row(s) exceed the float64 range on the response scale "
            f"and are returned as inf (first rows: {overflowed[:5].tolist()}). The location or scale prediction is "
            f"far outside the training range there.",
            UserWarning, stacklevel=3,
        )
    return response

The additive.py bridge (one-line delegates)

    def predict_quantiles(self, data: pd.DataFrame, taus: Sequence[float] = (0.05, 0.5, 0.95)) -> pd.DataFrame:
        return _distributional.predict_quantiles(self, data, taus)