# GammaLS


Gamma location-scale family for GAMLSS.


Usage

``` python
GammaLS()
```


[GammaLS](GammaLS.md#whittaker.GammaLS) extends the plain [Gamma](Gamma.md#whittaker.Gamma) family by letting both the mean `mu` and the coefficient of variation `sigma` vary smoothly with covariates, instead of assuming a fixed shape parameter. Use it for strictly positive, right-skewed responses where not only the typical magnitude but also the relative spread (coefficient of variation) changes systematically across the range of the predictors -- for example, cost or duration data whose relative volatility grows with the covariates rather than staying proportional to `mu` alone. Both parameters use the log link, keeping `mu > 0` and `sigma > 0`.


## Notes

[GammaLS](GammaLS.md#whittaker.GammaLS) parameterizes the Gamma distribution by its mean `mu > 0` and its coefficient of variation `sigma > 0`, where `sigma = 1 / sqrt(shape)` and `shape = alpha = 1 / sigma^2`. Both parameters use the log link:

 g\_{\mu}(\mu) = \log(\mu), \qquad g\_{\sigma}(\sigma) = \log(\sigma). 

The response density is the Gamma density with shape `alpha = 1/sigma^2` and rate `alpha/mu`:

 f(y \mid \mu, \sigma) = \frac{(\alpha/\mu)^{\alpha}}{\Gamma(\alpha)}\\ y^{\alpha - 1} \exp\\\left(-\frac{\alpha y}{\mu}\right), \qquad \alpha = \frac{1}{\sigma^{2}}. 

Because `sigma` is the coefficient of variation, \operatorname{Var}(Y) = \sigma^2 \mu^2, so this is a direct generalization of [Gamma](Gamma.md#whittaker.Gamma)'s variance function `V(mu) = mu^2` in which the proportionality constant `sigma^2` is itself allowed to depend on covariates.


## Examples

Fit a GAMLSS where both the mean and relative spread of a positive response vary smoothly:


``` python
import numpy as np
import whittaker as wk

rng = np.random.default_rng(0)
n = 300
x = np.linspace(0, 2 * np.pi, n)
mu = np.exp(1.0 + 0.4 * np.sin(x))
sigma = 0.2 + 0.15 * np.abs(np.cos(x))
shape = 1.0 / sigma**2
y = rng.gamma(shape, mu / shape)

data = {"x": x, "y": y}

model = wk.GAMLSS(
    formulas={"mu": "y ~ s(x)", "sigma": "y ~ s(x)"},
    family=wk.GammaLS(),
)
model.fit(data)
print(model.summary())
```


    GAMLSS fit summary
    ========================================
    Family: GammaLS(mu=log, sigma=log)
    N obs: 300
    Global deviance: 666.8366
    AIC: 690.3882
    BIC: 734.0031
    Log-likelihood: -333.4183
    Converged: True (5 iterations)

    --- mu ---
      EDF total: 5.93
      Smooth 1: edf = 4.93

    --- sigma ---
      EDF total: 5.84
      Smooth 1: edf = 4.84


## Attributes

| Name | Description |
|----|----|
| [parameter_names](#parameter_names) | Names of the distributional parameters modeled by [GammaLS](GammaLS.md#whittaker.GammaLS). |

------------------------------------------------------------------------


### parameter_names


Names of the distributional parameters modeled by [GammaLS](GammaLS.md#whittaker.GammaLS).


`parameter_names: tuple[str, …]`


## Methods

| Name | Description |
|----|----|
| [d2l_dtheta2()](#d2l_dtheta2) | Expected Fisher information for `mu` or `sigma`. |
| [dl_dtheta()](#dl_dtheta) | First derivative of the Gamma log-likelihood with respect to `mu` or `sigma`. |
| [initialize()](#initialize) | Starting values for `mu` and `sigma` from the raw response. |
| [link()](#link) | Apply the log link for `mu` or `sigma`. |
| [link_derivative()](#link_derivative) | Derivative of the log link for `mu` or `sigma`. |
| [link_inverse()](#link_inverse) | Apply the inverse log link for `mu` or `sigma`. |
| [log_likelihood()](#log_likelihood) | Total Gamma log-likelihood at the current `mu` and `sigma`. |
| [simulate()](#simulate) | Simulate response values from a Gamma distribution with the given `mu`, `sigma`. |

------------------------------------------------------------------------


### d2l_dtheta2()


Expected Fisher information for `mu` or `sigma`.


Usage

``` python
d2l_dtheta2(
    param,
    y,
    params,
)
```


With shape `alpha = 1/sigma^2` and trigamma function \psi_1:

 -E\\\left\[\frac{\partial^2 \ell}{\partial \mu^2}\right\] = \frac{\alpha}{\mu^2}, \qquad -E\\\left\[\frac{\partial^2 \ell}{\partial \sigma^2}\right\] = \frac{4}{\sigma^4}\big(\alpha\\\psi_1(\alpha) - 1\big). 

The `sigma` term is additionally floored at `_SIGMA_FLOOR` to keep the working weight strictly positive.


#### Parameters


`param: str`  
Either `"mu"` or `"sigma"`.

`y: NDArray`  
Observed response values, shape `(n,)` (unused, kept for interface consistency).

`params: dict[str, NDArray]`  
Current estimates with keys `"mu"` and `"sigma"`, each of shape `(n,)`.


#### Returns


`NDArray`  
Elementwise working weights, shape `(n,)`, guaranteed positive.


------------------------------------------------------------------------


### dl_dtheta()


First derivative of the Gamma log-likelihood with respect to `mu` or `sigma`.


Usage

``` python
dl_dtheta(
    param,
    y,
    params,
)
```


With shape `alpha = 1/sigma^2`, the score functions are

 \frac{\partial \ell}{\partial \mu} = \frac{\alpha}{\mu}\left(\frac{y}{\mu} - 1\right), \qquad \frac{\partial \ell}{\partial \sigma} = \frac{2}{\sigma^3}\left( \psi(\alpha) - \log\alpha - \log\frac{y}{\mu} - 1 + \frac{y}{\mu}\right), 

where \psi is the digamma function.


#### Parameters


`param: str`  
Either `"mu"` or `"sigma"`.

`y: NDArray`  
Observed response values, shape `(n,)`.

`params: dict[str, NDArray]`  
Current estimates with keys `"mu"` and `"sigma"`, each of shape `(n,)`.


#### Returns


`NDArray`  
Elementwise first derivatives, shape `(n,)`.


------------------------------------------------------------------------


### initialize()


Starting values for `mu` and `sigma` from the raw response.


Usage

``` python
initialize(y)
```


`mu` is initialized at the (floored) observed values themselves, and `sigma` at the sample coefficient of variation of `y` (clamped to at least `0.1`), constant across observations.


#### Parameters


`y: NDArray`  
Observed response values, shape `(n,)`.


#### Returns


`dict[str, NDArray]`  
`{"mu": y_safe.copy(), "sigma": <constant array of the sample CV>}`.


------------------------------------------------------------------------


### link()


Apply the log link for `mu` or `sigma`.


Usage

``` python
link(
    param,
    values,
)
```


Both parameters use the log link, so [values](TermsPredictionResult.md#whittaker.TermsPredictionResult.values) are first floored (`_MU_FLOOR` for `mu`, `_SIGMA_FLOOR` for `sigma`) to avoid `log(0)` before taking the logarithm.


#### Parameters


`param: str`  
Either `"mu"` or `"sigma"`.

`values: NDArray`  
Parameter values on the natural scale, shape `(n,)`.


#### Returns


`NDArray`  
`log(values)`, with [values](TermsPredictionResult.md#whittaker.TermsPredictionResult.values) clamped away from zero first.


------------------------------------------------------------------------


### link_derivative()


Derivative of the log link for `mu` or `sigma`.


Usage

``` python
link_derivative(
    param,
    values,
)
```


For the log link, d\eta/d\theta = 1/\theta for either parameter.


#### Parameters


`param: str`  
Either `"mu"` or `"sigma"`, used only to select the appropriate numerical floor.

`values: NDArray`  
Parameter values on the natural scale at which to evaluate the derivative, shape `(n,)`.


#### Returns


`NDArray`  
`1 / values`, with [values](TermsPredictionResult.md#whittaker.TermsPredictionResult.values) clamped away from zero first.


------------------------------------------------------------------------


### link_inverse()


Apply the inverse log link for `mu` or `sigma`.


Usage

``` python
link_inverse(
    param,
    eta,
)
```


`eta` is clipped to `[-_ETA_MAX, _ETA_MAX]` before exponentiating, guarding against overflow when the additive predictor drifts to extreme values during fitting.


#### Parameters


`param: str`  
Either `"mu"` or `"sigma"` (unused directly, since both share the log link).

`eta: NDArray`  
Additive predictor values, shape `(n,)`.


#### Returns


`NDArray`  
`exp(eta)` after clipping `eta` for numerical stability.


------------------------------------------------------------------------


### log_likelihood()


Total Gamma log-likelihood at the current `mu` and `sigma`.


Usage

``` python
log_likelihood(
    y,
    params,
)
```


Evaluated using the shape/rate parameterization `alpha = 1/sigma^2`, `rate = alpha/mu`, with `y` and `mu` floored at `_MU_FLOOR` and `sigma` floored at `_SIGMA_FLOOR` for numerical stability.


#### Parameters


`y: NDArray`  
Observed response values, shape `(n,)`.

`params: dict[str, NDArray]`  
Current estimates with keys `"mu"` and `"sigma"`, each of shape `(n,)`.


#### Returns


`float`  
The log-likelihood summed over all observations.


------------------------------------------------------------------------


### simulate()


Simulate response values from a Gamma distribution with the given `mu`, `sigma`.


Usage

``` python
simulate(
    params,
    rng,
)
```


Converts to the shape/scale parameterization (`shape = 1/sigma^2`, `scale = mu * sigma^2`) expected by `rng.gamma`.


#### Parameters


`params: dict[str, NDArray]`  
Current estimates with keys `"mu"` and `"sigma"`, each of shape `(n,)`.

`rng: np.random.Generator`  
A NumPy random generator used to draw the simulated values.


#### Returns


`NDArray`  
Simulated response values, shape `(n,)`.