# BetaLS


Beta family for GAMLSS with mean-precision parameterisation.


Usage

``` python
BetaLS()
```


[BetaLS](BetaLS.md#whittaker.BetaLS) extends the plain [Beta](Beta.md#whittaker.Beta) family by letting both the mean `mu` and the precision `phi` vary smoothly with covariates, rather than treating precision as a single estimated constant. Use it for responses strictly between 0 and 1 (rates, proportions, fractions) where not only the typical level but also how tightly the response clusters around that level changes across the range of the predictors -- for example, a proportion that becomes more variable in some regions of the covariate space and more tightly concentrated in others. The mean uses the logit link (as in [Beta](Beta.md#whittaker.Beta)) and the precision uses the log link, keeping `mu` in `(0, 1)` and `phi > 0`.


## Notes

The two parameters use distinct link functions:

 g\_{\mu}(\mu) = \log\\\left(\frac{\mu}{1-\mu}\right), \qquad g\_{\phi}(\phi) = \log(\phi). 

If `a = mu * phi` and `b = (1 - mu) * phi`, the response follows `y ~ Beta(a, b)` with density

 f(y \mid \mu, \phi) = \frac{y^{a-1}(1-y)^{b-1}}{B(a, b)}, \qquad a = \mu\phi,\\ \\ b = (1-\mu)\phi. 

As in [Beta](Beta.md#whittaker.Beta), larger `phi` concentrates the distribution more tightly around `mu` (`Var(Y) = mu(1-mu) / (1+phi)`), but here `phi` is itself modeled as a smooth function of covariates rather than a single scalar.


## Examples

Fit a GAMLSS where both the mean and precision of a proportion response vary smoothly:


``` python
import numpy as np
import whittaker as wk
from scipy.special import expit

rng = np.random.default_rng(0)
n = 300
x = np.linspace(0, 2 * np.pi, n)
mu = expit(np.sin(x))
phi = 10.0 + 15.0 * np.abs(np.cos(x))
y = rng.beta(mu * phi, (1 - mu) * phi)

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

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


    GAMLSS fit summary
    ========================================
    Family: BetaLS(mu=logit, phi=log)
    N obs: 300
    Global deviance: -548.1336
    AIC: -521.9058
    BIC: -473.3348
    Log-likelihood: 274.0668
    Converged: True (6 iterations)

    --- mu ---
      EDF total: 6.37
      Smooth 1: edf = 5.37

    --- phi ---
      EDF total: 6.74
      Smooth 1: edf = 5.74


## Attributes

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

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


### parameter_names


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


`parameter_names: tuple[str, …]`


## Methods

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

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


### d2l_dtheta2()


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


Usage

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


With `a = mu * phi`, `b = (1 - mu) * phi`, and trigamma function \psi_1:

 -E\\\left\[\frac{\partial^2 \ell}{\partial \mu^2}\right\] = \phi^2\big(\psi_1(a) + \psi_1(b)\big), \qquad -E\\\left\[\frac{\partial^2 \ell}{\partial \phi^2}\right\] = \mu^2 \psi_1(a) + (1-\mu)^2 \psi_1(b) - \psi_1(\phi). 


#### Parameters


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

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

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


#### Returns


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


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


### dl_dtheta()


First derivative of the Beta log-likelihood with respect to `mu` or `phi`.


Usage

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


With `a = mu * phi`, `b = (1 - mu) * phi`, `y* = logit(y)`, and `mu* = psi(a) - psi(b)` (where \psi is the digamma function), the score functions are

 \frac{\partial \ell}{\partial \mu} = \phi\\(y^{\*} - \mu^{\*}), \qquad \frac{\partial \ell}{\partial \phi} = \mu\\(y^{\*} - \mu^{\*}) + \psi(\phi) - \psi(b) + \log(1 - y). 


#### Parameters


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

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

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


#### Returns


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


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


### initialize()


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


Usage

``` python
initialize(y)
```


`mu` is initialized at the observed values themselves (clipped to `[0.01, 0.99]` to stay strictly inside the unit interval), and `phi` at a method-of-moments estimate of the precision derived from the sample mean and variance of `y` (clamped to at least `1.0`), constant across observations.


#### Parameters


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


#### Returns


`dict[str, NDArray]`  
`{"mu": y_safe.copy(), "phi": <constant array of the moment-based precision>}`.


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


### link()


Apply the link function for `mu` or `phi`.


Usage

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


`mu` uses the logit link, keeping its additive predictor unconstrained while `link_inverse` maps it back into `(0, 1)`; `phi` uses the log link, keeping it positive.


#### Parameters


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

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


#### Returns


`NDArray`  
`logit(values)` for `"mu"`, `log(values)` for `"phi"`.


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


### link_derivative()


Derivative of the link function for `mu` or `phi`.


Usage

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


For the logit link, d\eta/d\mu = 1/(\mu(1-\mu)); for the log link, d\eta/d\phi = 1/\phi.


#### Parameters


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

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


#### Returns


`NDArray`  
`1 / (values * (1 - values))` for `"mu"`, `1 / values` for `"phi"`.


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


### link_inverse()


Apply the inverse link for `mu` or `phi`.


Usage

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


Maps the additive predictor back to the natural scale: the logistic (`expit`) function for `mu`, restoring `(0, 1)`, and the exponential for `phi`, restoring positivity.


#### Parameters


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

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


#### Returns


`NDArray`  
`expit(eta)` for `"mu"`, `exp(eta)` for `"phi"`.


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


### log_likelihood()


Total Beta log-likelihood at the current `mu` and `phi`.


Usage

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


Evaluated using the shape parameterization `a = mu * phi`, `b = (1 - mu) * phi`.


#### Parameters


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

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


#### Returns


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


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


### simulate()


Simulate response values from a Beta distribution with the given `mu`, `phi`.


Usage

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


Converts to the shape parameterization (`a = mu * phi`, `b = (1 - mu) * phi`) expected by `rng.beta`.


#### Parameters


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

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


#### Returns


`NDArray`  
Simulated response values in `(0, 1)`, shape `(n,)`.