Fitting and Roadmap

roadmap

Development roadmap and future features for HeavyTails library.

This module contains placeholder functions and TODO items that will be automatically converted to GitHub Issues by the TODO workflow.

bootstrap_confidence_intervals

bootstrap_confidence_intervals(
    data,
    distribution,
    n_bootstrap=1000,
    confidence_level=0.95,
    seed=None,
)

Calculate bootstrap confidence intervals for distribution parameters.

Uses percentile bootstrap method to quantify uncertainty in MLE estimates.

Parameters:

Name	Type	Description	Default
data	list[float]	Sample data	required
distribution	str	Name of distribution to fit	required
n_bootstrap	int	Number of bootstrap samples (default: 1000)	1000
confidence_level	float	Confidence level, e.g., 0.95 for 95% CI (default: 0.95)	0.95
seed	int | None	Random seed for reproducibility (default: None)	None

Returns:

Type	Description
dict[str, tuple[float, float]]	Dictionary with parameter names as keys and (lower, upper) CI tuples as values

Examples:

>>> from heavytails import Pareto
>>> dist = Pareto(alpha=2.5, xm=1.0)
>>> data = dist.rvs(500, seed=42)
>>> ci = bootstrap_confidence_intervals(data, 'pareto', n_bootstrap=100, seed=42)
>>> 'alpha' in ci
True
>>> ci['alpha'][0] < 2.5 < ci['alpha'][1]  # Should contain true value
True

Source code in heavytails/roadmap.py

def bootstrap_confidence_intervals(
    data: list[float],
    distribution: str,
    n_bootstrap: int = 1000,
    confidence_level: float = 0.95,
    seed: int | None = None,
) -> dict[str, tuple[float, float]]:
    """
    Calculate bootstrap confidence intervals for distribution parameters.

    Uses percentile bootstrap method to quantify uncertainty in MLE estimates.

    Args:
        data: Sample data
        distribution: Name of distribution to fit
        n_bootstrap: Number of bootstrap samples (default: 1000)
        confidence_level: Confidence level, e.g., 0.95 for 95% CI (default: 0.95)
        seed: Random seed for reproducibility (default: None)

    Returns:
        Dictionary with parameter names as keys and (lower, upper) CI tuples as values

    Examples:
        >>> from heavytails import Pareto
        >>> dist = Pareto(alpha=2.5, xm=1.0)
        >>> data = dist.rvs(500, seed=42)
        >>> ci = bootstrap_confidence_intervals(data, 'pareto', n_bootstrap=100, seed=42)
        >>> 'alpha' in ci
        True
        >>> ci['alpha'][0] < 2.5 < ci['alpha'][1]  # Should contain true value
        True
    """
    if not data or len(data) == 0:
        raise ValueError("Data cannot be empty")

    if not 0 < confidence_level < 1:
        raise ValueError("Confidence level must be between 0 and 1")

    if n_bootstrap < 100:
        warnings.warn(
            "n_bootstrap < 100 may give unreliable confidence intervals", stacklevel=2
        )

    import random  # noqa: PLC0415

    # Set random seed for reproducibility
    if seed is not None:
        random.seed(seed)

    n = len(data)

    # Store bootstrap estimates
    bootstrap_estimates: dict[str, list[float]] = {
        param: [] for param in fit_mle(data, distribution)
    }

    # Perform bootstrap resampling
    for _ in range(n_bootstrap):
        # Resample with replacement
        bootstrap_sample = random.choices(data, k=n)

        try:
            # Fit distribution to bootstrap sample
            params = fit_mle(bootstrap_sample, distribution)

            # Store estimates
            for param_name, param_value in params.items():
                bootstrap_estimates[param_name].append(param_value)

        except (ValueError, RuntimeError):
            # Skip failed bootstrap samples
            continue

    # Calculate percentile confidence intervals
    alpha = 1 - confidence_level
    lower_percentile = 100 * (alpha / 2)
    upper_percentile = 100 * (1 - alpha / 2)

    confidence_intervals = {}

    for param_name, estimates in bootstrap_estimates.items():
        if len(estimates) < 10:
            warnings.warn(
                f"Too few successful bootstrap samples for {param_name}",
                stacklevel=2,
            )
            confidence_intervals[param_name] = (float("nan"), float("nan"))
            continue

        # Sort estimates
        sorted_estimates = sorted(estimates)

        # Calculate percentiles
        lower_idx = int(lower_percentile / 100 * len(sorted_estimates))
        upper_idx = int(upper_percentile / 100 * len(sorted_estimates))

        lower_bound = sorted_estimates[lower_idx]
        upper_bound = sorted_estimates[upper_idx]

        confidence_intervals[param_name] = (float(lower_bound), float(upper_bound))

    return confidence_intervals

fit_mle

fit_mle(data, distribution)

Fit distribution parameters using Maximum Likelihood Estimation.

Supports analytical and numerical MLE for all distributions in the library. For distributions without closed-form MLEs, uses scipy.optimize if available.

Parameters:

Name	Type	Description	Default
data	list[float]	Sample data to fit	required
distribution	str	Name of distribution (case-insensitive) Supported: 'pareto', 'lognormal', 'weibull', 'cauchy', 'studentt', 'exponential', 'frechet', 'generalizedpareto', 'burrxii', 'loglogistic', 'inversegamma', 'betaprime'	required

Returns:

Type	Description
dict[str, float]	Dictionary of fitted parameter names and values

Raises:

Type	Description
ValueError	If distribution is unknown or data is invalid
ImportError	If scipy is required but not available

Examples:

>>> import random
>>> from heavytails import Pareto
>>> dist = Pareto(alpha=2.5, xm=1.0)
>>> data = dist.rvs(1000, seed=42)
>>> params = fit_mle(data, 'pareto')
>>> abs(params['alpha'] - 2.5) < 0.2  # Should be close
True

Source code in heavytails/roadmap.py

def fit_mle(data: list[float], distribution: str) -> dict[str, float]:
    """
    Fit distribution parameters using Maximum Likelihood Estimation.

    Supports analytical and numerical MLE for all distributions in the library.
    For distributions without closed-form MLEs, uses scipy.optimize if available.

    Args:
        data: Sample data to fit
        distribution: Name of distribution (case-insensitive)
            Supported: 'pareto', 'lognormal', 'weibull', 'cauchy', 'studentt',
            'exponential', 'frechet', 'generalizedpareto', 'burrxii',
            'loglogistic', 'inversegamma', 'betaprime'

    Returns:
        Dictionary of fitted parameter names and values

    Raises:
        ValueError: If distribution is unknown or data is invalid
        ImportError: If scipy is required but not available

    Examples:
        >>> import random
        >>> from heavytails import Pareto
        >>> dist = Pareto(alpha=2.5, xm=1.0)
        >>> data = dist.rvs(1000, seed=42)
        >>> params = fit_mle(data, 'pareto')
        >>> abs(params['alpha'] - 2.5) < 0.2  # Should be close
        True
    """
    if not data or len(data) == 0:
        raise ValueError("Data cannot be empty")

    if any(not math.isfinite(x) for x in data):
        raise ValueError("Data contains non-finite values")

    dist_lower = distribution.lower()

    # Dictionary of MLE estimators
    estimators = {
        "pareto": _fit_pareto_mle,
        "lognormal": _fit_lognormal_mle,
        "weibull": _fit_weibull_mle,
        "cauchy": _fit_cauchy_mle,
        "studentt": _fit_studentt_mle,
        "exponential": _fit_exponential_mle,
        "frechet": _fit_frechet_mle,
        "generalizedpareto": _fit_gpd_mle,
        "burrxii": _fit_burrxii_mle,
        "loglogistic": _fit_loglogistic_mle,
        "inversegamma": _fit_inversegamma_mle,
        "betaprime": _fit_betaprime_mle,
    }

    if dist_lower not in estimators:
        available = ", ".join(sorted(estimators.keys()))
        raise ValueError(
            f"MLE not implemented for '{distribution}'. Available: {available}"
        )

    return estimators[dist_lower](data)

model_comparison

model_comparison(data, distributions)

Compare distribution fits using information criteria.

Computes AIC and BIC for each distribution and ranks them. Lower values indicate better fit (penalized by model complexity).

Parameters:

Name	Type	Description	Default
data	list[float]	Sample data	required
distributions	list[str]	List of distribution names to compare	required

Returns:

Type	Description
dict[str, dict[str, Any]]	Dictionary with results for each distribution containing: - params: Fitted parameters - log_likelihood: Log-likelihood value - AIC: Akaike Information Criterion - BIC: Bayesian Information Criterion - rank_AIC: Rank by AIC (1 = best) - rank_BIC: Rank by BIC (1 = best)

Examples:

>>> from heavytails import Pareto
>>> dist = Pareto(alpha=2.5, xm=1.0)
>>> data = dist.rvs(1000, seed=42)
>>> results = model_comparison(data, ['pareto', 'lognormal', 'weibull'])
>>> results['pareto']['rank_AIC']  # Pareto should rank best
1

Source code in heavytails/roadmap.py

def model_comparison(
    data: list[float], distributions: list[str]
) -> dict[str, dict[str, Any]]:
    """
    Compare distribution fits using information criteria.

    Computes AIC and BIC for each distribution and ranks them.
    Lower values indicate better fit (penalized by model complexity).

    Args:
        data: Sample data
        distributions: List of distribution names to compare

    Returns:
        Dictionary with results for each distribution containing:
            - params: Fitted parameters
            - log_likelihood: Log-likelihood value
            - AIC: Akaike Information Criterion
            - BIC: Bayesian Information Criterion
            - rank_AIC: Rank by AIC (1 = best)
            - rank_BIC: Rank by BIC (1 = best)

    Examples:
        >>> from heavytails import Pareto
        >>> dist = Pareto(alpha=2.5, xm=1.0)
        >>> data = dist.rvs(1000, seed=42)
        >>> results = model_comparison(data, ['pareto', 'lognormal', 'weibull'])
        >>> results['pareto']['rank_AIC']  # Pareto should rank best
        1
    """
    if not data or len(data) == 0:
        raise ValueError("Data cannot be empty")

    n = len(data)
    results = {}

    for dist_name in distributions:
        try:
            # Fit distribution
            params = fit_mle(data, dist_name)

            # Calculate log-likelihood
            log_likelihood = _calculate_log_likelihood(data, dist_name, params)

            # Number of parameters
            k = len(params)

            # Calculate information criteria
            aic = 2 * k - 2 * log_likelihood
            bic = k * math.log(n) - 2 * log_likelihood

            results[dist_name] = {
                "params": params,
                "log_likelihood": log_likelihood,
                "AIC": aic,
                "BIC": bic,
                "n_params": k,
            }

        except (ValueError, ImportError) as e:
            warnings.warn(
                f"Failed to fit {dist_name}: {e}",
                stacklevel=2,
            )
            results[dist_name] = {
                "params": None,
                "log_likelihood": float("-inf"),
                "AIC": float("inf"),
                "BIC": float("inf"),
                "n_params": 0,
                "error": str(e),
            }

    # Rank models by AIC and BIC
    valid_results = {
        k: v for k, v in results.items() if v["log_likelihood"] != float("-inf")
    }

    if valid_results:
        aic_sorted = sorted(valid_results.items(), key=lambda x: float(x[1]["AIC"]))  # type: ignore[arg-type]
        bic_sorted = sorted(valid_results.items(), key=lambda x: float(x[1]["BIC"]))  # type: ignore[arg-type]

        for rank, (dist_name, _) in enumerate(aic_sorted, 1):
            results[dist_name]["rank_AIC"] = rank

        for rank, (dist_name, _) in enumerate(bic_sorted, 1):
            results[dist_name]["rank_BIC"] = rank

    return results

robust_hill_estimator

robust_hill_estimator(data, k=None, bias_correction=True)

Improved Hill estimator with bias correction and stability checks.

Implements bias-corrected Hill estimator with automatic k selection and diagnostic information for assessing estimate reliability.

Parameters:

Name	Type	Description	Default
data	list[float]	Sample data (should be heavy-tailed)	required
k	int | None	Number of top order statistics to use. If None, automatically selected.	None
bias_correction	bool	Apply second-order bias correction (default: True)	True

Returns:

Type	Description
dict[str, float | int | bool | str]	Dictionary containing: - gamma: Tail index estimate (gamma = 1/alpha for Pareto) - alpha: Shape parameter estimate (alpha = 1/gamma) - k_used: Number of order statistics used - bias_corrected: Whether bias correction was applied - n: Sample size - reliability: Quality indicator ('good', 'fair', 'poor')

Examples:

>>> from heavytails import Pareto
>>> dist = Pareto(alpha=2.5, xm=1.0)
>>> data = dist.rvs(1000, seed=42)
>>> result = robust_hill_estimator(data)
>>> abs(result['alpha'] - 2.5) < 0.5  # Should be close
True
>>> result['reliability'] in ['good', 'fair', 'poor']
True

References

Dekkers, A. L., Einmahl, J. H., & De Haan, L. (1989). A moment estimator for the index of an extreme-value distribution. Annals of Statistics, 17(4), 1833-1855.

Source code in heavytails/roadmap.py

def robust_hill_estimator(
    data: list[float], k: int | None = None, bias_correction: bool = True
) -> dict[str, float | int | bool | str]:
    """
    Improved Hill estimator with bias correction and stability checks.

    Implements bias-corrected Hill estimator with automatic k selection
    and diagnostic information for assessing estimate reliability.

    Args:
        data: Sample data (should be heavy-tailed)
        k: Number of top order statistics to use. If None, automatically selected.
        bias_correction: Apply second-order bias correction (default: True)

    Returns:
        Dictionary containing:
            - gamma: Tail index estimate (gamma = 1/alpha for Pareto)
            - alpha: Shape parameter estimate (alpha = 1/gamma)
            - k_used: Number of order statistics used
            - bias_corrected: Whether bias correction was applied
            - n: Sample size
            - reliability: Quality indicator ('good', 'fair', 'poor')

    Examples:
        >>> from heavytails import Pareto
        >>> dist = Pareto(alpha=2.5, xm=1.0)
        >>> data = dist.rvs(1000, seed=42)
        >>> result = robust_hill_estimator(data)
        >>> abs(result['alpha'] - 2.5) < 0.5  # Should be close
        True
        >>> result['reliability'] in ['good', 'fair', 'poor']
        True

    References:
        Dekkers, A. L., Einmahl, J. H., & De Haan, L. (1989).
        A moment estimator for the index of an extreme-value distribution.
        Annals of Statistics, 17(4), 1833-1855.
    """
    n = len(data)

    # Sample size check
    if n < 50:
        raise ValueError("Sample size too small (n < 50) for reliable Hill estimation")

    if n < 200:
        warnings.warn(
            "Hill estimator may be unreliable for n < 200. "
            "Consider collecting more data.",
            stacklevel=2,
        )

    # Automatic k selection if not provided
    if k is None:
        k = _select_optimal_k(data)

    # Validate k
    if not (5 < k < n // 2):
        k = max(5, min(k, n // 2 - 1))
        warnings.warn(f"k adjusted to valid range: k = {k}", stacklevel=2)

    # Basic Hill estimate
    gamma_basic = hill_estimator(data, k)

    if not bias_correction:
        alpha = 1.0 / gamma_basic if gamma_basic > 0 else float("inf")
        return {
            "gamma": gamma_basic,
            "alpha": alpha,
            "k_used": k,
            "bias_corrected": False,
            "n": n,
            "reliability": _assess_reliability(n, k),
        }

    # Apply bias correction (Dekkers-Einmahl-de Haan)
    sorted_data = sorted(data, reverse=True)
    x_k = sorted_data[k]

    # Calculate higher-order moments for bias correction
    logs = [math.log(sorted_data[i] / x_k) for i in range(k)]
    M1 = sum(logs) / k
    M2 = sum(log_val**2 for log_val in logs) / k

    # Bias-corrected estimate
    if M2 > M1**2:
        rho = 1.0 - 0.5 * (1.0 - M1**2 / M2) ** -1
        bias = rho * gamma_basic / (1.0 - rho)

        # Apply correction with safeguards
        gamma_corrected = gamma_basic - bias

        # Ensure corrected estimate is reasonable
        if gamma_corrected <= 0 or gamma_corrected > 2 * gamma_basic:
            warnings.warn(
                "Bias correction yielded unreasonable estimate, using basic Hill",
                stacklevel=2,
            )
            gamma_final = gamma_basic
        else:
            gamma_final = gamma_corrected
    else:
        # Cannot apply bias correction
        warnings.warn("Insufficient moment variation for bias correction", stacklevel=2)
        gamma_final = gamma_basic

    alpha = 1.0 / gamma_final if gamma_final > 0 else float("inf")

    return {
        "gamma": gamma_final,
        "alpha": alpha,
        "k_used": k,
        "bias_corrected": True,
        "n": n,
        "reliability": _assess_reliability(n, k),
    }