insurance_distributional_glm.families

146class GamlssFamily(ABC):
147    """
148    Abstract base for GAMLSS distribution families.
149
150    Each subclass implements a specific distribution and its relationship
151    to the RS fitting algorithm. The key quantities are:
152      - log_likelihood: total log-likelihood (sum over obs)
153      - dl_deta: d(log L)/d(eta_k) — the score for the IRLS working response
154      - d2l_deta2: E[d^2(log L)/d(eta_k)^2] — negative for IRLS weights
155
156    These drive the weighted least squares step in the RS algorithm.
157    """
158
159    @property
160    @abstractmethod
161    def param_names(self) -> List[str]:
162        """Ordered list of distribution parameter names, e.g. ['mu', 'sigma']."""
163
164    @property
165    @abstractmethod
166    def default_links(self) -> Dict[str, Link]:
167        """Default link function for each parameter."""
168
169    @abstractmethod
170    def log_likelihood(self, y: np.ndarray, params: Dict[str, np.ndarray]) -> np.ndarray:
171        """
172        Per-observation log-likelihood.
173
174        Returns array of shape (n,), one value per observation.
175        The fitting engine sums these to get total log-likelihood.
176        """
177
178    @abstractmethod
179    def dl_deta(
180        self,
181        y: np.ndarray,
182        params: Dict[str, np.ndarray],
183        param_name: str,
184    ) -> np.ndarray:
185        """
186        d(log L_i)/d(eta_k) for parameter param_name.
187
188        This is the derivative of log-likelihood with respect to the linear
189        predictor eta_k (not theta_k). Used to construct the working response
190        in the IRLS step.
191
192        By chain rule: d(log L)/d(eta_k) = d(log L)/d(theta_k) * d(theta_k)/d(eta_k)
193
194        Raises ValueError for unknown param_name.
195        """
196
197    @abstractmethod
198    def d2l_deta2(
199        self,
200        y: np.ndarray,
201        params: Dict[str, np.ndarray],
202        param_name: str,
203    ) -> np.ndarray:
204        """
205        Expected -d^2(log L_i)/d(eta_k)^2 for parameter param_name.
206
207        Returns positive values (this is the negative expected second derivative).
208        Used as IRLS weights. If analytic expectation is hard, use the observed
209        Fisher information (negative of actual second derivative).
210
211        Raises ValueError for unknown param_name.
212        """
213
214    @abstractmethod
215    def starting_values(self, y: np.ndarray) -> Dict[str, np.ndarray]:
216        """
217        Initial parameter estimates from the marginal distribution of y.
218
219        Called once before the RS loop begins. Should be robust — moment
220        estimates are fine here; precision doesn't matter much.
221        """
222
223    def _check_param(self, param_name: str) -> None:
224        """Raise ValueError if param_name is not a known parameter."""
225        if param_name not in self.param_names:
226            raise ValueError(
227                f"Unknown parameter '{param_name}' for {self.__class__.__name__}. "
228                f"Known parameters: {self.param_names}"
229            )
230
231    def validate_params(self, params: Dict[str, np.ndarray]) -> None:
232        """Optional validation hook. Raise ValueError if params are out of range."""
233        pass
234
235    def __repr__(self) -> str:
236        return f"{self.__class__.__name__}(links={{{', '.join(f'{k}: {v.name}' for k, v in self.default_links.items())}}})"

26class Link(ABC):
27    """Abstract link function g: parameter space -> real line."""
28
29    @abstractmethod
30    def link(self, theta: np.ndarray) -> np.ndarray:
31        """Forward: eta = g(theta)."""
32
33    @abstractmethod
34    def inverse(self, eta: np.ndarray) -> np.ndarray:
35        """Inverse: theta = g^{-1}(eta)."""
36
37    @abstractmethod
38    def deriv(self, theta: np.ndarray) -> np.ndarray:
39        """Derivative: d(eta)/d(theta) = g'(theta)."""
40
41    def inverse_deriv(self, eta: np.ndarray) -> np.ndarray:
42        """d(theta)/d(eta) = 1 / g'(g^{-1}(eta))."""
43        theta = self.inverse(eta)
44        return 1.0 / self.deriv(theta)
45
46    @property
47    @abstractmethod
48    def name(self) -> str:
49        pass

52class LogLink(Link):
53    """Log link: eta = log(theta), theta = exp(eta). Used for positive parameters."""
54
55    def link(self, theta: np.ndarray) -> np.ndarray:
56        return np.log(theta)
57
58    def inverse(self, eta: np.ndarray) -> np.ndarray:
59        return np.exp(np.clip(eta, -30, 30))
60
61    def deriv(self, theta: np.ndarray) -> np.ndarray:
62        return 1.0 / theta
63
64    def inverse_deriv(self, eta: np.ndarray) -> np.ndarray:
65        return np.exp(np.clip(eta, -30, 30))
66
67    @property
68    def name(self) -> str:
69        return "log"

72class LogitLink(Link):
73    """Logit link: eta = log(theta/(1-theta)), theta = expit(eta). Used for probabilities."""
74
75    def link(self, theta: np.ndarray) -> np.ndarray:
76        theta_clipped = np.clip(theta, 1e-10, 1 - 1e-10)
77        return np.log(theta_clipped / (1.0 - theta_clipped))
78
79    def inverse(self, eta: np.ndarray) -> np.ndarray:
80        return 1.0 / (1.0 + np.exp(-np.clip(eta, -30, 30)))
81
82    def deriv(self, theta: np.ndarray) -> np.ndarray:
83        theta_clipped = np.clip(theta, 1e-10, 1 - 1e-10)
84        return 1.0 / (theta_clipped * (1.0 - theta_clipped))
85
86    def inverse_deriv(self, eta: np.ndarray) -> np.ndarray:
87        p = self.inverse(eta)
88        return p * (1.0 - p)
89
90    @property
91    def name(self) -> str:
92        return "logit"

 95class IdentityLink(Link):
 96    """Identity link: eta = theta. Used for location parameters on the original scale."""
 97
 98    def link(self, theta: np.ndarray) -> np.ndarray:
 99        return theta.copy()
100
101    def inverse(self, eta: np.ndarray) -> np.ndarray:
102        return eta.copy()
103
104    def deriv(self, theta: np.ndarray) -> np.ndarray:
105        return np.ones_like(theta)
106
107    def inverse_deriv(self, eta: np.ndarray) -> np.ndarray:
108        return np.ones_like(eta)
109
110    @property
111    def name(self) -> str:
112        return "identity"

115class SqrtLink(Link):
116    """Square root link: eta = sqrt(theta), theta = eta^2."""
117
118    def link(self, theta: np.ndarray) -> np.ndarray:
119        return np.sqrt(np.clip(theta, 0, None))
120
121    def inverse(self, eta: np.ndarray) -> np.ndarray:
122        return np.clip(eta, 0, None) ** 2
123
124    def deriv(self, theta: np.ndarray) -> np.ndarray:
125        return 0.5 / np.sqrt(np.clip(theta, 1e-10, None))
126
127    def inverse_deriv(self, eta: np.ndarray) -> np.ndarray:
128        return 2.0 * np.clip(eta, 0, None)
129
130    @property
131    def name(self) -> str:
132        return "sqrt"

 28class Gamma(GamlssFamily):
 29    """
 30    Gamma distribution with mean mu and coefficient of variation sigma.
 31
 32    PDF: f(y) = (1/sigma^2)^(1/sigma^2) * y^(1/sigma^2 - 1) * exp(-y/(mu*sigma^2))
 33                / (Gamma(1/sigma^2) * (mu*sigma^2)^(1/sigma^2))
 34
 35    mu    > 0 : mean (log link)
 36    sigma > 0 : coefficient of variation, sd/mean (log link)
 37
 38    Insurance use: claim severity. sigma models heterogeneity of variance
 39    — more risky segments have higher CV.
 40    """
 41
 42    @property
 43    def param_names(self):
 44        return ["mu", "sigma"]
 45
 46    @property
 47    def default_links(self):
 48        return {"mu": LogLink(), "sigma": LogLink()}
 49
 50    def log_likelihood(self, y: np.ndarray, params: Dict[str, np.ndarray]) -> np.ndarray:
 51        mu = params["mu"]
 52        sigma = params["sigma"]
 53        nu = 1.0 / sigma**2  # shape parameter
 54        return (
 55            nu * np.log(nu)
 56            - gammaln(nu)
 57            + (nu - 1.0) * np.log(y)
 58            - nu * np.log(mu)
 59            - nu * y / mu
 60        )
 61
 62    def dl_deta(self, y, params, param_name):
 63        self._check_param(param_name)
 64        mu = params["mu"]
 65        sigma = params["sigma"]
 66        nu = 1.0 / sigma**2
 67
 68        if param_name == "mu":
 69            link = self.default_links["mu"]
 70            dl_dtheta = nu * (y - mu) / mu**2
 71            dtheta_deta = link.inverse_deriv(link.link(mu))
 72            return dl_dtheta * dtheta_deta
 73
 74        else:  # sigma
 75            from scipy.special import digamma
 76            link = self.default_links["sigma"]
 77            dl_dnu = np.log(nu) + 1.0 - digamma(nu) + np.log(y / mu) - y / mu
 78            dl_dsigma = dl_dnu * (-2.0 / sigma**3)
 79            dsigma_deta = link.inverse_deriv(link.link(sigma))
 80            return dl_dsigma * dsigma_deta
 81
 82    def d2l_deta2(self, y, params, param_name):
 83        self._check_param(param_name)
 84        mu = params["mu"]
 85        sigma = params["sigma"]
 86        nu = 1.0 / sigma**2
 87
 88        if param_name == "mu":
 89            link = self.default_links["mu"]
 90            d2l_dmu2 = nu / mu**2
 91            dtheta_deta = link.inverse_deriv(link.link(mu))
 92            return d2l_dmu2 * dtheta_deta**2
 93
 94        else:  # sigma
 95            from scipy.special import polygamma
 96            link = self.default_links["sigma"]
 97            d2l_dnu2 = polygamma(1, nu) - 1.0 / nu
 98            dnu_dsigma = -2.0 / sigma**3
 99            d2l_dsigma2 = d2l_dnu2 * dnu_dsigma**2
100            dsigma_deta = link.inverse_deriv(link.link(sigma))
101            return d2l_dsigma2 * dsigma_deta**2
102
103    def starting_values(self, y):
104        mu_init = np.full(len(y), np.mean(y))
105        sigma_init = np.full(len(y), np.std(y) / np.mean(y))
106        sigma_init = np.clip(sigma_init, 0.01, 5.0)
107        return {"mu": mu_init, "sigma": sigma_init}

110class LogNormal(GamlssFamily):
111    """
112    Log-Normal distribution with log-scale mean mu and log-scale SD sigma.
113
114    If Y ~ LogNormal(mu, sigma), then log(Y) ~ Normal(mu, sigma).
115    Note: mu here is the mean of log(Y), not E[Y].
116    E[Y] = exp(mu + sigma^2/2).
117
118    mu    : real (identity link — this is already on log scale)
119    sigma > 0 : log-scale standard deviation (log link)
120
121    Insurance use: severity when log(claims) is approximately normal.
122    Heavier right tail than Gamma for the same mean.
123    """
124
125    @property
126    def param_names(self):
127        return ["mu", "sigma"]
128
129    @property
130    def default_links(self):
131        return {"mu": IdentityLink(), "sigma": LogLink()}
132
133    def log_likelihood(self, y, params):
134        mu = params["mu"]
135        sigma = params["sigma"]
136        log_y = np.log(np.clip(y, 1e-300, None))
137        return (
138            -np.log(sigma)
139            - 0.5 * np.log(2.0 * np.pi)
140            - log_y
141            - 0.5 * ((log_y - mu) / sigma) ** 2
142        )
143
144    def dl_deta(self, y, params, param_name):
145        self._check_param(param_name)
146        mu = params["mu"]
147        sigma = params["sigma"]
148        log_y = np.log(np.clip(y, 1e-300, None))
149        resid = (log_y - mu) / sigma
150
151        if param_name == "mu":
152            link = self.default_links["mu"]
153            dl_dmu = resid / sigma
154            dmu_deta = link.inverse_deriv(link.link(mu))
155            return dl_dmu * dmu_deta
156
157        else:  # sigma
158            link = self.default_links["sigma"]
159            dl_dsigma = (resid**2 - 1.0) / sigma
160            dsigma_deta = link.inverse_deriv(link.link(sigma))
161            return dl_dsigma * dsigma_deta
162
163    def d2l_deta2(self, y, params, param_name):
164        self._check_param(param_name)
165        mu = params["mu"]
166        sigma = params["sigma"]
167
168        if param_name == "mu":
169            link = self.default_links["mu"]
170            d2l_dmu2 = 1.0 / sigma**2
171            dmu_deta = link.inverse_deriv(link.link(mu))
172            return d2l_dmu2 * dmu_deta**2
173
174        else:  # sigma
175            link = self.default_links["sigma"]
176            d2l_dsigma2 = 2.0 / sigma**2
177            dsigma_deta = link.inverse_deriv(link.link(sigma))
178            return d2l_dsigma2 * dsigma_deta**2
179
180    def starting_values(self, y):
181        log_y = np.log(np.clip(y, 1e-300, None))
182        mu_init = np.full(len(y), np.mean(log_y))
183        sigma_init = np.full(len(y), max(np.std(log_y), 0.01))
184        return {"mu": mu_init, "sigma": sigma_init}

187class InverseGaussian(GamlssFamily):
188    """
189    Inverse Gaussian distribution with mean mu and dispersion sigma.
190
191    Standard parameterisation: if lambda = 1/sigma^2 is the shape (precision)
192    parameter, then:
193      PDF: f(y) = sqrt(lambda/(2*pi*y^3)) * exp(-lambda*(y-mu)^2/(2*mu^2*y))
194      Var[Y] = mu^3 / lambda = sigma^2 * mu^3
195
196    mu    > 0 : mean (log link)
197    sigma > 0 : dispersion, Var = sigma^2 * mu^3 (log link)
198
199    Log-likelihood:
200      ll = -log(sigma) - 0.5*log(2*pi) - 1.5*log(y) - (y-mu)^2/(2*sigma^2*mu^2*y)
201
202    Insurance use: liability severity with heavy right tails. The variance
203    grows faster with mu than Gamma (mu^3 vs mu^2), suitable for large losses.
204    """
205
206    @property
207    def param_names(self):
208        return ["mu", "sigma"]
209
210    @property
211    def default_links(self):
212        return {"mu": LogLink(), "sigma": LogLink()}
213
214    def log_likelihood(self, y, params):
215        mu = params["mu"]
216        sigma = params["sigma"]
217        y_safe = np.clip(y, 1e-300, None)
218        return (
219            -np.log(sigma)
220            - 0.5 * np.log(2.0 * np.pi)
221            - 1.5 * np.log(y_safe)
222            - (y - mu) ** 2 / (2.0 * sigma**2 * mu**2 * y_safe)
223        )
224
225    def dl_deta(self, y, params, param_name):
226        self._check_param(param_name)
227        mu = params["mu"]
228        sigma = params["sigma"]
229        y_safe = np.clip(y, 1e-300, None)
230
231        if param_name == "mu":
232            link = self.default_links["mu"]
233            # d(ll)/d(mu):
234            # ll contains -(y-mu)^2/(2*sigma^2*mu^2*y)
235            # Expand: -(y^2/mu^2 - 2y/mu + 1)/(2*sigma^2*y)
236            # d/dmu = (2y^2/mu^3 - 2y/mu^2)/(2*sigma^2*y) = (y-mu)/(sigma^2*mu^3)
237            dl_dmu = (y - mu) / (sigma**2 * mu**3)
238            dmu_deta = link.inverse_deriv(link.link(mu))
239            return dl_dmu * dmu_deta
240
241        else:  # sigma
242            link = self.default_links["sigma"]
243            # d(ll)/d(sigma) = -1/sigma + (y-mu)^2/(sigma^3*mu^2*y)
244            dl_dsigma = -1.0 / sigma + (y - mu)**2 / (sigma**3 * mu**2 * y_safe)
245            dsigma_deta = link.inverse_deriv(link.link(sigma))
246            return dl_dsigma * dsigma_deta
247
248    def d2l_deta2(self, y, params, param_name):
249        self._check_param(param_name)
250        mu = params["mu"]
251        sigma = params["sigma"]
252
253        if param_name == "mu":
254            link = self.default_links["mu"]
255            # E[-d^2(ll)/dmu^2] = 1/(sigma^2*mu^3)
256            d2l_dmu2 = 1.0 / (sigma**2 * mu**3)
257            dmu_deta = link.inverse_deriv(link.link(mu))
258            return d2l_dmu2 * dmu_deta**2
259
260        else:  # sigma
261            link = self.default_links["sigma"]
262            # E[-d^2 ll/d sigma^2] = 2/sigma^2
263            d2l_dsigma2 = 2.0 / sigma**2
264            dsigma_deta = link.inverse_deriv(link.link(sigma))
265            return d2l_dsigma2 * dsigma_deta**2
266
267    def starting_values(self, y):
268        mu_init = np.full(len(y), np.mean(y))
269        sigma_sq = max(np.var(y) / np.mean(y)**3, 1e-6)
270        sigma_init = np.full(len(y), np.sqrt(sigma_sq))
271        return {"mu": mu_init, "sigma": sigma_init}