gamma

gamma ¶

Gamma distribution for ngboost-lightning.

Gamma ¶

Gamma(params: NDArray[floating])

Bases: Distribution

Gamma distribution with log-link parameterization.

Internal parameters are [log_alpha, log_beta] where alpha = exp(log_alpha) is the shape parameter and beta = exp(log_beta) is the rate parameter. The log-links keep both parameters positive during unconstrained boosting.

Note on Fisher Information

For Gamma(log_alpha, log_beta), the Fisher Information is non-diagonal. The off-diagonal entries are -alpha. This distribution relies on the base class np.linalg.solve for the natural gradient — no fast-path override.

ATTRIBUTE	DESCRIPTION
`n_params`	Always 2 (log_alpha and log_beta).
`alpha`	Shape parameter values, shape `[n_samples]`. TYPE: `NDArray[floating]`
`beta`	Rate parameter values, shape `[n_samples]`. TYPE: `NDArray[floating]`

Construct Gamma from internal parameters.

PARAMETER	DESCRIPTION
`params`	Internal parameters, shape `[n_samples, 2]`. Column 0 is log(alpha), column 1 is log(beta). TYPE: `NDArray[floating]`

Source code in ngboost_lightning/distributions/gamma.py

def __init__(self, params: NDArray[np.floating]) -> None:
    """Construct Gamma from internal parameters.

    Args:
        params: Internal parameters, shape ``[n_samples, 2]``.
            Column 0 is log(alpha), column 1 is log(beta).
    """
    self.log_alpha: NDArray[np.floating] = params[:, 0]
    self.log_beta: NDArray[np.floating] = params[:, 1]
    self.alpha: NDArray[np.floating] = np.exp(self.log_alpha)
    self.beta: NDArray[np.floating] = np.exp(self.log_beta)
    self._dist = sp_gamma(a=self.alpha, scale=1.0 / self.beta)
    self._params = params

fit `staticmethod` ¶

fit(
    y: NDArray[floating],
    sample_weight: NDArray[floating] | None = None,
) -> NDArray[floating]

Estimate initial (log_alpha, log_beta) via method of moments.

PARAMETER	DESCRIPTION
`y`	Target values, shape `[n_samples]`. Must be positive. TYPE: `NDArray[floating]`
`sample_weight`	Per-sample weights, shape `[n_samples]`. TYPE: `NDArray[floating] \| None` DEFAULT: `None`

RETURNS	DESCRIPTION
`NDArray[floating]`	Parameter vector `[log(alpha), log(beta)]`, shape `[2]`.

Source code in ngboost_lightning/distributions/gamma.py

@staticmethod
def fit(
    y: NDArray[np.floating],
    sample_weight: NDArray[np.floating] | None = None,
) -> NDArray[np.floating]:
    """Estimate initial (log_alpha, log_beta) via method of moments.

    Args:
        y: Target values, shape ``[n_samples]``. Must be positive.
        sample_weight: Per-sample weights, shape ``[n_samples]``.

    Returns:
        Parameter vector ``[log(alpha), log(beta)]``, shape ``[2]``.
    """
    mean_y = float(np.average(y, weights=sample_weight))
    var_y = float(np.average((y - mean_y) ** 2, weights=sample_weight))
    beta = np.maximum(mean_y / max(var_y, 1e-10), 1e-6)
    alpha = np.maximum(mean_y * beta, 1e-6)
    return np.array([np.log(alpha), np.log(beta)])

score ¶

score(y: NDArray[floating]) -> NDArray[floating]

Per-sample negative log-likelihood.

PARAMETER	DESCRIPTION
`y`	Observed target values, shape `[n_samples]`. TYPE: `NDArray[floating]`

RETURNS	DESCRIPTION
`NDArray[floating]`	NLL values, shape `[n_samples]`.

Source code in ngboost_lightning/distributions/gamma.py

def score(self, y: NDArray[np.floating]) -> NDArray[np.floating]:
    """Per-sample negative log-likelihood.

    Args:
        y: Observed target values, shape ``[n_samples]``.

    Returns:
        NLL values, shape ``[n_samples]``.
    """
    return -self._dist.logpdf(y)

d_score ¶

d_score(y: NDArray[floating]) -> NDArray[floating]

Analytical gradient of NLL w.r.t. [log_alpha, log_beta].

Derivation

NLL = -(alpha-1)log(y) + betay - alpha*log(beta) + gammaln(alpha)

d(NLL)/d(alpha) = -log(y) - log(beta) + digamma(alpha) d(NLL)/d(log_alpha) = d(NLL)/d(alpha) * alpha = alpha * (digamma(alpha) - log(beta) - log(y))

d(NLL)/d(beta) = y - alpha/beta d(NLL)/d(log_beta) = d(NLL)/d(beta) * beta = beta*y - alpha

PARAMETER	DESCRIPTION
`y`	Observed target values, shape `[n_samples]`. TYPE: `NDArray[floating]`

RETURNS	DESCRIPTION
`NDArray[floating]`	Gradient array, shape `[n_samples, 2]`.

Source code in ngboost_lightning/distributions/gamma.py

def d_score(self, y: NDArray[np.floating]) -> NDArray[np.floating]:
    """Analytical gradient of NLL w.r.t. [log_alpha, log_beta].

    Derivation:
        NLL = -(alpha-1)*log(y) + beta*y - alpha*log(beta) + gammaln(alpha)

        d(NLL)/d(alpha) = -log(y) - log(beta) + digamma(alpha)
        d(NLL)/d(log_alpha) = d(NLL)/d(alpha) * alpha
            = alpha * (digamma(alpha) - log(beta) - log(y))

        d(NLL)/d(beta) = y - alpha/beta
        d(NLL)/d(log_beta) = d(NLL)/d(beta) * beta
            = beta*y - alpha

    Args:
        y: Observed target values, shape ``[n_samples]``.

    Returns:
        Gradient array, shape ``[n_samples, 2]``.
    """
    n = len(y)
    grad = np.empty((n, 2))
    grad[:, 0] = self.alpha * (digamma(self.alpha) - self.log_beta - np.log(y))
    grad[:, 1] = self.beta * y - self.alpha
    return grad

metric ¶

metric() -> NDArray[floating]

Fisher Information for (log_alpha, log_beta) parameterization.

The standard-parameterization Fisher Information for Gamma(alpha, beta) is::

FI = [[trigamma(alpha), -1 / beta], [-1 / beta, alpha / beta ^ 2]]

Applying the Jacobian J = diag(alpha, beta) for the log-link::

FI_log = J ^ T @ FI @ J

gives::

FI_log = [[alpha ^ 2 * trigamma(alpha), -alpha], [-alpha, alpha]]

RETURNS	DESCRIPTION
`NDArray[floating]`	FI tensor, shape `[n_samples, 2, 2]`.

Source code in ngboost_lightning/distributions/gamma.py

def metric(self) -> NDArray[np.floating]:
    """Fisher Information for (log_alpha, log_beta) parameterization.

    The standard-parameterization Fisher Information for Gamma(alpha, beta)
    is::

        FI = [[trigamma(alpha), -1 / beta], [-1 / beta, alpha / beta ^ 2]]

    Applying the Jacobian ``J = diag(alpha, beta)`` for the log-link::

        FI_log = J ^ T @ FI @ J

    gives::

        FI_log = [[alpha ^ 2 * trigamma(alpha), -alpha], [-alpha, alpha]]

    Returns:
        FI tensor, shape ``[n_samples, 2, 2]``.
    """
    n = len(self.alpha)
    fi = np.empty((n, 2, 2))
    trigamma = polygamma(1, self.alpha)
    fi[:, 0, 0] = self.alpha**2 * trigamma
    fi[:, 0, 1] = -self.alpha
    fi[:, 1, 0] = -self.alpha
    fi[:, 1, 1] = self.alpha
    return fi

mean ¶

mean() -> NDArray[floating]

Conditional mean (point prediction): alpha / beta.

RETURNS	DESCRIPTION
`NDArray[floating]`	Mean values, shape `[n_samples]`.

Source code in ngboost_lightning/distributions/gamma.py

def mean(self) -> NDArray[np.floating]:
    """Conditional mean (point prediction): alpha / beta.

    Returns:
        Mean values, shape ``[n_samples]``.
    """
    return self.alpha / self.beta

sample ¶

sample(n: int) -> NDArray[floating]

Draw n samples per distribution instance.

PARAMETER	DESCRIPTION
`n`	Number of samples to draw. TYPE: `int`

RETURNS	DESCRIPTION
`NDArray[floating]`	Samples, shape `[n, n_samples]`.

Source code in ngboost_lightning/distributions/gamma.py

def sample(self, n: int) -> NDArray[np.floating]:
    """Draw n samples per distribution instance.

    Args:
        n: Number of samples to draw.

    Returns:
        Samples, shape ``[n, n_samples]``.
    """
    return self._dist.rvs(size=(n, len(self)))

cdf ¶

cdf(y: NDArray[floating]) -> NDArray[floating]

Cumulative distribution function.

PARAMETER	DESCRIPTION
`y`	Values at which to evaluate the CDF. TYPE: `NDArray[floating]`

RETURNS	DESCRIPTION
`NDArray[floating]`	CDF values, same shape as `y`.

Source code in ngboost_lightning/distributions/gamma.py

def cdf(self, y: NDArray[np.floating]) -> NDArray[np.floating]:
    """Cumulative distribution function.

    Args:
        y: Values at which to evaluate the CDF.

    Returns:
        CDF values, same shape as ``y``.
    """
    return self._dist.cdf(y)

ppf ¶

ppf(q: NDArray[floating]) -> NDArray[floating]

Percent point function (inverse CDF / quantile function).

PARAMETER	DESCRIPTION
`q`	Quantiles, values in [0, 1]. TYPE: `NDArray[floating]`

RETURNS	DESCRIPTION
`NDArray[floating]`	Values at the given quantiles, same shape as `q`.

Source code in ngboost_lightning/distributions/gamma.py

def ppf(self, q: NDArray[np.floating]) -> NDArray[np.floating]:
    """Percent point function (inverse CDF / quantile function).

    Args:
        q: Quantiles, values in [0, 1].

    Returns:
        Values at the given quantiles, same shape as ``q``.
    """
    return self._dist.ppf(q)

logpdf ¶

logpdf(y: NDArray[floating]) -> NDArray[floating]

Log probability density function.

PARAMETER	DESCRIPTION
`y`	Values at which to evaluate. TYPE: `NDArray[floating]`

RETURNS	DESCRIPTION
`NDArray[floating]`	Log-density values, same shape as `y`.

Source code in ngboost_lightning/distributions/gamma.py

def logpdf(self, y: NDArray[np.floating]) -> NDArray[np.floating]:
    """Log probability density function.

    Args:
        y: Values at which to evaluate.

    Returns:
        Log-density values, same shape as ``y``.
    """
    return self._dist.logpdf(y)

crps_score ¶

crps_score(y: NDArray[floating]) -> NDArray[floating]

Per-sample CRPS for Gamma.

Closed form (Scheuerer & Hamill 2015): CRPS = y*(2*F(y;a,b) - 1) - (a/b)*(2*F(y;a+1,b) - 1) - 1/(b * Beta(a, 0.5))

where a = alpha (shape), b = beta (rate), F(y;a,b) is the Gamma CDF with shape a and rate b.

PARAMETER	DESCRIPTION
`y`	Observed target values, shape `[n_samples]`. TYPE: `NDArray[floating]`

RETURNS	DESCRIPTION
`NDArray[floating]`	CRPS values, shape `[n_samples]`.

Source code in ngboost_lightning/distributions/gamma.py

def crps_score(self, y: NDArray[np.floating]) -> NDArray[np.floating]:
    """Per-sample CRPS for Gamma.

    Closed form (Scheuerer & Hamill 2015):
        ``CRPS = y*(2*F(y;a,b) - 1)
                 - (a/b)*(2*F(y;a+1,b) - 1)
                 - 1/(b * Beta(a, 0.5))``

    where ``a = alpha`` (shape), ``b = beta`` (rate),
    ``F(y;a,b)`` is the Gamma CDF with shape ``a`` and rate ``b``.

    Args:
        y: Observed target values, shape ``[n_samples]``.

    Returns:
        CRPS values, shape ``[n_samples]``.
    """
    a = self.alpha
    b = self.beta
    # F(y; shape=a, rate=b) = F(y; shape=a, scale=1/b)
    cdf_a = sp_gamma.cdf(y, a=a, scale=1.0 / b)
    cdf_a1 = sp_gamma.cdf(y, a=a + 1.0, scale=1.0 / b)
    return (
        y * (2.0 * cdf_a - 1.0)
        - (a / b) * (2.0 * cdf_a1 - 1.0)
        - 1.0 / (b * sp_beta(a, 0.5))
    )

crps_d_score ¶

crps_d_score(y: NDArray[floating]) -> NDArray[floating]

Gradient of CRPS w.r.t. [log_alpha, log_beta].

Uses central finite differences on the internal parameters for robustness. The analytical gradient of the Gamma CDF w.r.t. the shape parameter involves the derivative of the regularized incomplete gamma function, which is numerically delicate. Finite differences are accurate to O(eps^2) and fast enough for the boosting loop since CRPS evaluation is cheap.

PARAMETER	DESCRIPTION
`y`	Observed target values, shape `[n_samples]`. TYPE: `NDArray[floating]`

RETURNS	DESCRIPTION
`NDArray[floating]`	Gradient array, shape `[n_samples, 2]`.

Source code in ngboost_lightning/distributions/gamma.py

def crps_d_score(self, y: NDArray[np.floating]) -> NDArray[np.floating]:
    """Gradient of CRPS w.r.t. [log_alpha, log_beta].

    Uses central finite differences on the internal parameters for
    robustness.  The analytical gradient of the Gamma CDF w.r.t. the
    shape parameter involves the derivative of the regularized
    incomplete gamma function, which is numerically delicate.  Finite
    differences are accurate to ``O(eps^2)`` and fast enough for
    the boosting loop since CRPS evaluation is cheap.

    Args:
        y: Observed target values, shape ``[n_samples]``.

    Returns:
        Gradient array, shape ``[n_samples, 2]``.
    """
    eps = 1e-5
    n = len(y)
    grad = np.empty((n, 2))
    params = self._params.copy()

    for k in range(2):
        params_plus = params.copy()
        params_plus[:, k] += eps
        params_minus = params.copy()
        params_minus[:, k] -= eps
        score_plus = type(self)(params_plus).crps_score(y)
        score_minus = type(self)(params_minus).crps_score(y)
        grad[:, k] = (score_plus - score_minus) / (2.0 * eps)

    return grad

crps_metric ¶

crps_metric() -> NDArray[floating]

Riemannian metric for CRPS natural gradient (MC estimate).

Uses a Monte Carlo estimate: sample from the distribution, compute CRPS gradients, and average the outer product.

RETURNS	DESCRIPTION
`NDArray[floating]`	Metric tensor, shape `[n_samples, 2, 2]`.

Source code in ngboost_lightning/distributions/gamma.py

def crps_metric(self) -> NDArray[np.floating]:
    """Riemannian metric for CRPS natural gradient (MC estimate).

    Uses a Monte Carlo estimate: sample from the distribution, compute
    CRPS gradients, and average the outer product.

    Returns:
        Metric tensor, shape ``[n_samples, 2, 2]``.
    """
    n_mc = 50
    n = len(self.alpha)
    met = np.zeros((n, 2, 2))
    rng = np.random.default_rng(42)
    for _ in range(n_mc):
        y_sample = rng.gamma(self.alpha, 1.0 / self.beta)
        g = self.crps_d_score(y_sample)
        # Outer product: g[:, :, None] @ g[:, None, :]
        met += np.einsum("ij,ik->ijk", g, g)
    met /= n_mc
    return met

gamma

gamma ¶

Gamma ¶

fit staticmethod ¶

score ¶

d_score ¶

metric ¶

mean ¶

sample ¶

cdf ¶

ppf ¶

logpdf ¶

crps_score ¶

crps_d_score ¶

crps_metric ¶

fit `staticmethod` ¶