API Reference¶
This page documents the public API and internal engines in magli.
Public Unified API¶
The top-level GAM (aliased as LinGAM) and GAMCore classes dynamically dispatch
to either the Wood (backend='wood', default) or Eilers-Marx (backend='eilers') engine.
Backend-specific classes are also directly importable from magli as
WoodLinGAM, WoodGAMCore, EilersLinGAM, and EilersGAMCore.
Eilers & Marx Engine (magli.eilers)¶
B-spline additive models with discrete GCV grid search and soft shape constraints following Eilers & Marx (1996).
Classes¶
Modules¶
magli.eilers._bspline¶
Vectorized De Boor B-spline basis construction.
- magli.eilers._bspline.b_spline_basis(x_in: ArrayLike, n_splines: int, spline_order: int, edge_knots: ndarray, periodic: bool = False) ndarray[source]¶
Construct B-spline basis matrix using vectorized De Boor recursion.
- Parameters:
x_in (array-like, shape (n,)) – Predictor values at which to evaluate the basis.
n_splines (int) – Number of basis functions K. For periodic splines this is the full count before wrapping; the returned matrix has
K - spline_ordercolumns.spline_order (int) – Order k of the B-spline (polynomial degree).
edge_knots (np.ndarray, shape (2,)) – Domain boundaries [min, max] for mapping data to [0, 1].
periodic (bool, default False) – If True, wrap the first
spline_orderbasis functions to enforce periodic continuity across f(0) = f(1).
- Returns:
bases – B-spline basis matrix with linear extrapolation beyond edge knots.
- Return type:
np.ndarray, shape (n, K) or (n, K_eff) for periodic
1"""Vectorized De Boor B-spline basis construction."""
2
3import numpy as np
4from numpy.typing import ArrayLike
5
6
7def b_spline_basis(
8 x_in: ArrayLike,
9 n_splines: int,
10 spline_order: int,
11 edge_knots: np.ndarray,
12 periodic: bool = False,
13) -> np.ndarray:
14 """Construct B-spline basis matrix using vectorized De Boor recursion.
15
16 Parameters
17 ----------
18 x_in : array-like, shape (n,)
19 Predictor values at which to evaluate the basis.
20 n_splines : int
21 Number of basis functions K. For periodic splines this is the
22 *full* count before wrapping; the returned matrix has
23 ``K - spline_order`` columns.
24 spline_order : int
25 Order k of the B-spline (polynomial degree).
26 edge_knots : np.ndarray, shape (2,)
27 Domain boundaries [min, max] for mapping data to [0, 1].
28 periodic : bool, default False
29 If True, wrap the first ``spline_order`` basis functions
30 to enforce periodic continuity across f(0) = f(1).
31
32 Returns
33 -------
34 bases : np.ndarray, shape (n, K) or (n, K_eff) for periodic
35 B-spline basis matrix with linear extrapolation beyond edge knots.
36 """
37 K = n_splines
38 k = spline_order
39
40 if K <= k:
41 raise ValueError(
42 f"n_splines ({K}) must be greater than spline_order ({k}). "
43 f"Use at least {k + 1}."
44 )
45
46 offset = edge_knots[0]
47 scale = edge_knots[1] - edge_knots[0]
48 if scale == 0.0:
49 scale = 1.0
50
51 boundary_knots = np.linspace(0.0, 1.0, 1 + K - k)
52 diff = boundary_knots[1] - boundary_knots[0]
53
54 x = (np.asarray(x_in).ravel() - offset) / scale
55 x = np.r_[x, 0.0, 1.0]
56
57 if periodic:
58 x = x % 1.0
59 x[-2] = 0.0
60 x[-1] = 1.0
61
62 x_extrap_l = x < 0.0
63 x_extrap_r = x > 1.0
64 x_interp = ~(x_extrap_l | x_extrap_r)
65 x = np.atleast_2d(x).T
66
67 aug = np.arange(1, k + 1) * diff
68 aug_knots = np.r_[-aug[::-1], boundary_knots, 1.0 + aug]
69 aug_knots[-1] += 1e-9
70
71 bases = ((x >= aug_knots[:-1]) & (x < aug_knots[1:])).astype(float)
72 bases[-1] = bases[-2][::-1]
73
74 n_active = len(aug_knots) - 1
75 prev_bases = None
76
77 for m in range(2, k + 2):
78 n_active -= 1
79 num_l = (x - aug_knots[:n_active]) * bases[:, :n_active]
80 denom_l = aug_knots[m - 1:n_active + m - 1] - aug_knots[:n_active]
81 left = num_l / denom_l
82
83 num_r = (aug_knots[m:n_active + m] - x) * bases[:, 1:n_active + 1]
84 denom_r = aug_knots[m:n_active + m] - aug_knots[1:n_active + 1]
85 right = num_r / denom_r
86
87 prev_bases = bases[-2:].copy()
88 bases = left + right
89
90 if (np.any(x_extrap_l) or np.any(x_extrap_r)) and k > 0 and not periodic:
91 bases[~x_interp] = 0.0
92 denom_l = aug_knots[k:-1] - aug_knots[:-k - 1]
93 left_g = prev_bases[:, :-1] / denom_l
94 denom_r = aug_knots[k + 1:] - aug_knots[1:-k]
95 right_g = prev_bases[:, 1:] / denom_r
96 grads = k * (left_g - right_g)
97
98 if np.any(x_extrap_l):
99 bases[x_extrap_l] = grads[0] * x[x_extrap_l] + bases[-2]
100 if np.any(x_extrap_r):
101 bases[x_extrap_r] = grads[1] * (x[x_extrap_r] - 1.0) + bases[-1]
102
103 result = bases[:-2]
104
105 if periodic:
106 n_wrap = k
107 if n_wrap > 0 and n_wrap < K:
108 K_eff = K - n_wrap
109 wrapped = result[:, :K_eff].copy()
110 wrapped[:, :n_wrap] += result[:, K_eff:K_eff + n_wrap]
111 return wrapped
112
113 return result
magli.eilers._penalty¶
Penalty matrix builders for GAM regularisation.
- magli.eilers._penalty.second_diff_penalty(k: int, lam: float) ndarray[source]¶
Second-difference penalty D2^T D2 scaled by lam.
- magli.eilers._penalty.first_diff_penalty(k: int, lam: float) ndarray[source]¶
First-difference penalty D1^T D1 scaled by lam.
- magli.eilers._penalty.circular_second_diff_penalty(k: int, lam: float) ndarray[source]¶
Circular second-difference penalty for periodic splines.
- magli.eilers._penalty.ridge_penalty(k: int, lam: float) ndarray[source]¶
L2 ridge penalty lam * I_k.
- magli.eilers._penalty.shape_constraint_penalty(coef: ndarray, constraint: str, lam_con: float = 1000000000.0, l2_con: float = 0.001) ndarray[source]¶
Asymmetric active-set penalty for monotonicity and convexity constraints.
1"""Penalty matrix builders for GAM regularisation."""
2
3import numpy as np
4
5
6def second_diff_penalty(k: int, lam: float) -> np.ndarray:
7 """Second-difference penalty D2^T D2 scaled by lam."""
8 D2 = np.diff(np.eye(k), n=2, axis=0)
9 return lam * (D2.T @ D2)
10
11
12def first_diff_penalty(k: int, lam: float) -> np.ndarray:
13 """First-difference penalty D1^T D1 scaled by lam."""
14 D1 = np.diff(np.eye(k), n=1, axis=0)
15 return lam * (D1.T @ D1)
16
17
18def circular_second_diff_penalty(k: int, lam: float) -> np.ndarray:
19 """Circular second-difference penalty for periodic splines."""
20 I = np.eye(k)
21 D2 = np.roll(I, -1, axis=1) - 2.0 * I + np.roll(I, 1, axis=1)
22 return lam * (D2.T @ D2)
23
24
25def ridge_penalty(k: int, lam: float) -> np.ndarray:
26 """L2 ridge penalty lam * I_k."""
27 return lam * np.eye(k)
28
29
30def shape_constraint_penalty(
31 coef: np.ndarray,
32 constraint: str,
33 lam_con: float = 1e9,
34 l2_con: float = 1e-3,
35) -> np.ndarray:
36 """Asymmetric active-set penalty for monotonicity and convexity constraints."""
37 k = len(coef)
38 if k <= 1:
39 return np.zeros((k, k))
40 if constraint in ("mono_inc", "mono_dec"):
41 D = np.diff(np.eye(k), n=1, axis=0)
42 d1 = np.diff(coef)
43 mask = (d1 < 0) if constraint == "mono_inc" else (d1 > 0)
44 elif constraint in ("convex", "concave"):
45 if k <= 2:
46 return np.zeros((k, k))
47 D = np.diff(np.eye(k), n=2, axis=0)
48 d2 = np.diff(coef, n=2)
49 mask = (d2 < 0) if constraint == "convex" else (d2 > 0)
50 else:
51 return np.zeros((k, k))
52
53 if not np.any(mask):
54 return np.zeros((k, k))
55 Dm = D[mask]
56 C = lam_con * (Dm.T @ Dm)
57 C.flat[:: k + 1] += l2_con
58 return C
59
magli.eilers._solver¶
Penalized least-squares solver via SVD data augmentation.
- magli.eilers._solver.solve_pirls(X: ndarray, y: ndarray, P: ndarray) Tuple[ndarray, ndarray, ndarray][source]¶
Solve penalized least squares: (X^T X + P) beta = X^T y.
1"""Penalized least-squares solver via SVD data augmentation."""
2
3from typing import Tuple
4
5import numpy as np
6
7
8def solve_pirls(
9 X: np.ndarray,
10 y: np.ndarray,
11 P: np.ndarray,
12) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
13 """Solve penalized least squares: (X^T X + P) beta = X^T y."""
14 n, m = X.shape
15 eps = np.finfo(float).eps
16 E = np.linalg.cholesky(np.diag(np.ones(m) * np.sqrt(eps)) + P).T
17 Q, R = np.linalg.qr(X)
18 min_nm = min(m, n)
19 U, d, Vt = np.linalg.svd(np.vstack([R, E]), full_matrices=False)
20 U1 = U[:min_nm, :]
21 B_solve = (Vt.T / d) @ U1.T @ Q.T
22 return B_solve @ y, U1, B_solve
magli.eilers._terms¶
GAM term classes: spline, tensor, factor, and linear.
- magli.eilers._terms.instantiate_terms(configs: List[Dict], spline_order: int) List[_Term][source]¶
Create term objects from parsed config dicts.
1"""GAM term classes: spline, tensor, factor, and linear."""
2
3from typing import Dict, List, Optional, Tuple, Union
4
5import numpy as np
6
7from .._base import _khatri_rao, _kron_sum
8from ._bspline import b_spline_basis
9from ._penalty import (
10 circular_second_diff_penalty,
11 ridge_penalty,
12 second_diff_penalty,
13 shape_constraint_penalty,
14)
15
16_CONSTRAINTS = frozenset({
17 None, 'mono_inc', 'mono_dec', 'convex', 'concave', 'periodic',
18})
19_SHAPE_CONSTRAINTS = frozenset({'mono_inc', 'mono_dec', 'convex', 'concave'})
20
21
22class _Term:
23 """Base class for all GAM terms."""
24
25 n_coefs: int = 0
26 sp: Optional[List[float]] = None
27
28 def compile(self, x: np.ndarray) -> None:
29 """One-time setup using training data."""
30
31 def build_penalties(self) -> List[np.ndarray]:
32 return []
33
34
35class _SplineTerm(_Term):
36 """Univariate penalised B-spline smooth term."""
37
38 __slots__ = (
39 'feature', 'n_splines', 'lam', 'spline_order',
40 '_edge_knots', 'constraint', 'n_coefs',
41 )
42
43 def __init__(
44 self,
45 feature: int,
46 n_splines: int = 10,
47 lam: float = 1.0,
48 spline_order: int = 3,
49 constraint: Optional[str] = None,
50 ) -> None:
51 if constraint not in _CONSTRAINTS:
52 raise ValueError(
53 f"Unknown constraint '{constraint}'. "
54 f"Choose from: {sorted(c for c in _CONSTRAINTS if c is not None)}."
55 )
56 min_k = (2 * spline_order) if constraint == 'periodic' else (spline_order + 1)
57 if n_splines < min_k:
58 raise ValueError(
59 f"n_splines ({n_splines}) must be at least {min_k} "
60 f"for spline_order ({spline_order})"
61 + (" with periodic constraint." if constraint == 'periodic' else ".")
62 )
63 self.feature = feature
64 self.n_splines = n_splines
65 self.lam = lam
66 self.spline_order = spline_order
67 self.constraint = constraint
68 self.n_coefs = (
69 n_splines - spline_order
70 if constraint == 'periodic'
71 else n_splines
72 )
73 self._edge_knots: Optional[np.ndarray] = None
74
75 def compile(self, x: np.ndarray) -> None:
76 col = np.asarray(x[:, self.feature]).ravel()
77 self._edge_knots = np.array([col.min(), col.max()])
78
79 def build_columns(self, x: np.ndarray) -> np.ndarray:
80 return b_spline_basis(
81 x[:, self.feature],
82 self.n_splines,
83 self.spline_order,
84 self._edge_knots, # type: ignore[arg-type]
85 periodic=self.constraint == 'periodic',
86 )
87
88 def build_penalty(self, lam: Optional[float] = None) -> np.ndarray:
89 l = lam if lam is not None else self.lam
90 K_eff = self.n_coefs
91 if self.constraint == 'periodic':
92 return circular_second_diff_penalty(K_eff, l)
93 return second_diff_penalty(K_eff, l)
94
95 @property
96 def has_shape_constraint(self) -> bool:
97 return self.constraint in _SHAPE_CONSTRAINTS
98
99 def build_constraint_penalty(self, coef_sub: np.ndarray) -> np.ndarray:
100 if self.constraint in _SHAPE_CONSTRAINTS:
101 return shape_constraint_penalty(coef_sub, self.constraint)
102 return np.zeros((self.n_coefs, self.n_coefs))
103
104 def penalty_bases(self) -> List[Tuple[np.ndarray, int]]:
105 return [(self.build_penalty(1.0), 1)]
106
107
108class _TensorTerm(_Term):
109 """Tensor-product B-spline interaction term (2+D)."""
110
111 __slots__ = (
112 'features', 'd', 'n_splines_per', 'lam_per', 'lam',
113 'spline_order', '_edge_knots', 'n_coefs',
114 )
115
116 def __init__(
117 self,
118 features: List[int],
119 n_splines: Union[int, List[int]] = 5,
120 lam: Union[float, List[float]] = 1.0,
121 spline_order: int = 3,
122 ) -> None:
123 self.features = list(features)
124 self.d = len(self.features)
125 if self.d < 2:
126 raise ValueError("Tensor term requires at least 2 features.")
127
128 if isinstance(n_splines, int):
129 self.n_splines_per = [n_splines] * self.d
130 elif len(n_splines) != self.d:
131 raise ValueError(
132 f"n_splines length ({len(n_splines)}) must match "
133 f"number of features ({self.d})."
134 )
135 else:
136 self.n_splines_per = list(n_splines)
137
138 if isinstance(lam, (int, float)):
139 self.lam_per = [float(lam)] * self.d
140 elif len(lam) != self.d:
141 raise ValueError(
142 f"lam length ({len(lam)}) must match "
143 f"number of features ({self.d})."
144 )
145 else:
146 self.lam_per = [float(v) for v in lam]
147
148 self.lam = self.lam_per[0]
149 self.spline_order = spline_order
150 self.n_coefs = int(np.prod(self.n_splines_per))
151 self._edge_knots: List[np.ndarray] = []
152
153 def compile(self, x: np.ndarray) -> None:
154 self._edge_knots = [
155 np.array([np.min(x[:, f]), np.max(x[:, f])]) for f in self.features
156 ]
157
158 def build_columns(self, x: np.ndarray) -> np.ndarray:
159 return _khatri_rao([
160 b_spline_basis(
161 x[:, feat],
162 self.n_splines_per[i],
163 self.spline_order,
164 self._edge_knots[i],
165 )
166 for i, feat in enumerate(self.features)
167 ])
168
169 def build_penalty(self, lam: Optional[List[float]] = None) -> np.ndarray:
170 lams = lam if lam is not None else self.lam_per
171 if isinstance(lams, (int, float)):
172 lams = [float(lams)] * self.d
173 return _kron_sum([
174 second_diff_penalty(self.n_splines_per[i], lams[i])
175 for i in range(self.d)
176 ])
177
178 def penalty_bases(self) -> List[Tuple[np.ndarray, int]]:
179 sizes = list(self.n_splines_per)
180 bases: List[Tuple[np.ndarray, int]] = []
181 for i, sz in enumerate(sizes):
182 S_i = second_diff_penalty(sz, 1.0)
183 I_pre = np.eye(int(np.prod(sizes[:i])))
184 I_post = np.eye(int(np.prod(sizes[i + 1:])))
185 bases.append((np.kron(np.kron(I_pre, S_i), I_post), 1))
186 return bases
187
188
189class _FactorTerm(_Term):
190 """Categorical / factor term with one-hot or dummy encoding."""
191
192 __slots__ = ('feature', 'lam', 'coding', '_levels', 'n_coefs')
193
194 def __init__(
195 self,
196 feature: int,
197 lam: float = 1.0,
198 coding: str = 'one-hot',
199 ) -> None:
200 self.feature = int(feature)
201 self.lam = lam
202 self.coding = coding
203 self._levels: Optional[np.ndarray] = None
204 self.n_coefs = 0
205
206 def compile(self, x: np.ndarray) -> None:
207 self._levels = np.unique(np.asarray(x[:, self.feature]).ravel())
208 self.n_coefs = (
209 len(self._levels)
210 if self.coding == 'one-hot'
211 else len(self._levels) - 1
212 )
213 if self.n_coefs <= 0:
214 raise ValueError(f"Factor term feature {self.feature} has no levels.")
215
216 def build_columns(self, x: np.ndarray) -> np.ndarray:
217 if self._levels is None:
218 self.compile(x)
219 col = np.asarray(x[:, self.feature]).ravel()
220 B = (col[:, None] == self._levels[None, :]).astype(float)
221 return B[:, 1:] if self.coding == 'dummy' else B
222
223 def build_penalty(self, lam: Optional[float] = None) -> np.ndarray:
224 return ridge_penalty(self.n_coefs, lam if lam is not None else self.lam)
225
226 def penalty_bases(self) -> List[Tuple[np.ndarray, int]]:
227 return [(np.eye(self.n_coefs), 1)]
228
229
230class _LinearTerm(_Term):
231 """Linear (unpenalised by default) term."""
232
233 __slots__ = ('feature', 'lam', 'n_coefs')
234
235 def __init__(self, feature: int, lam: float = 0.0) -> None:
236 self.feature = int(feature)
237 self.lam = lam
238 self.n_coefs = 1
239
240 def build_columns(self, x: np.ndarray) -> np.ndarray:
241 return np.asarray(x[:, self.feature], dtype=float).ravel().reshape(-1, 1)
242
243 def build_penalty(self, lam: Optional[float] = None) -> np.ndarray:
244 return np.array([[lam if lam is not None else self.lam]], dtype=float)
245
246 def penalty_bases(self) -> List[Tuple[np.ndarray, int]]:
247 return [(np.array([[1.0]]), 1)]
248
249
250def instantiate_terms(
251 configs: List[Dict],
252 spline_order: int,
253) -> List[_Term]:
254 """Create term objects from parsed config dicts."""
255 terms: List[_Term] = []
256 for cfg in configs:
257 t = cfg['type']
258 feats = cfg['features']
259 kw = dict(cfg['kwargs'])
260 term_order = kw.pop('spline_order', spline_order)
261 k_val = kw.pop('n_splines', kw.pop('k', None))
262 lam_val = kw.pop('lam', kw.pop('sp', None))
263
264 if t == 's':
265 terms.append(_SplineTerm(
266 feats[0],
267 k_val if k_val is not None else 10,
268 lam_val if lam_val is not None else 1.0,
269 term_order,
270 constraint=kw.pop('constraint', None),
271 ))
272 elif t == 'te':
273 terms.append(_TensorTerm(
274 feats,
275 k_val if k_val is not None else 5,
276 lam_val if lam_val is not None else 1.0,
277 term_order,
278 ))
279 elif t == 'f':
280 terms.append(_FactorTerm(
281 feats[0],
282 lam_val if lam_val is not None else 1.0,
283 kw.pop('coding', 'one-hot'),
284 ))
285 elif t == 'l':
286 terms.append(_LinearTerm(
287 feats[0],
288 lam_val if lam_val is not None else 0.0,
289 ))
290 else:
291 raise ValueError(f"Unknown term type: {t}")
292
293 if kw:
294 raise ValueError(
295 f"Unrecognised kwargs for {t}(): {list(kw.keys())}"
296 )
297 return terms
magli.eilers._lingam¶
Single-term penalised B-spline GAM (Eilers & Marx backend).
- class magli.eilers._lingam.LinGAM(n_splines: int = 10, lam: float = 1.0, spline_order: int = 3, constraint: str | None = None)[source]¶
Penalised B-spline GAM for a single smooth term y = f(x) + eps.
1"""Single-term penalised B-spline GAM (Eilers & Marx backend)."""
2
3from typing import Optional, Tuple
4
5import numpy as np
6from numpy.typing import ArrayLike
7
8from ._gam import GAMCore
9from ._terms import _SplineTerm
10
11
12class LinGAM(GAMCore):
13 """Penalised B-spline GAM for a single smooth term y = f(x) + eps."""
14
15 def __init__(
16 self,
17 n_splines: int = 10,
18 lam: float = 1.0,
19 spline_order: int = 3,
20 constraint: Optional[str] = None,
21 ) -> None:
22 if lam < 0:
23 raise ValueError("Smoothing parameter 'lam' must be non-negative.")
24 if spline_order < 1:
25 raise ValueError("spline_order must be >= 1.")
26 super().__init__(
27 f"s(0, n_splines={int(n_splines)}, lam={float(lam)}, constraint={constraint!r})",
28 spline_order=int(spline_order),
29 fit_intercept=True,
30 )
31
32 @property
33 def _term(self) -> _SplineTerm:
34 return self._terms[0] # type: ignore[return-value]
35
36 @property
37 def n_splines(self) -> int:
38 return self._term.n_splines
39
40 @n_splines.setter
41 def n_splines(self, val: int) -> None:
42 val = int(val)
43 self._term.n_splines = val
44 self._term.n_coefs = (
45 val - self.spline_order
46 if self.constraint == "periodic"
47 else val
48 )
49 self._term_configs[0]["kwargs"]["n_splines"] = val
50
51 @property
52 def lam(self) -> float:
53 return self._term.lam
54
55 @lam.setter
56 def lam(self, val: float) -> None:
57 self._term.lam = float(val)
58 self._term_configs[0]["kwargs"]["lam"] = float(val)
59
60 @property
61 def constraint(self) -> Optional[str]:
62 return self._term.constraint
63
64 @property
65 def n_coefs(self) -> int:
66 """Effective number of spline coefficients (reduced for periodic)."""
67 return self._term.n_coefs
68
69 @property
70 def edge_knots_(self) -> Optional[np.ndarray]:
71 return self._term._edge_knots
72
73 @property
74 def design_matrix_(self) -> Optional[np.ndarray]:
75 return self._design_matrix_
76
77 @property
78 def knots_(self) -> Optional[np.ndarray]:
79 ek = self.edge_knots_
80 if ek is None:
81 return None
82 scale = (ek[1] - ek[0]) or 1.0
83 return (
84 np.linspace(0.0, 1.0, 1 + self.n_splines - self.spline_order)
85 * scale
86 + ek[0]
87 )
88
89 def __repr__(self) -> str:
90 status = "fitted" if self.is_fitted else "unfitted"
91 cons = f", constraint='{self.constraint}'" if self.constraint else ""
92 return (
93 f"LinGAM(n_splines={self.n_splines}, lam={self.lam:.4g}, "
94 f"spline_order={self.spline_order}{cons}) [{status}]"
95 )
96
97 def _diagnostics_title(self) -> str:
98 cons = f", constraint='{self.constraint}'" if self.constraint else ""
99 return f"LinGAM Diagnostics (K={self.n_splines}, lam={self.lam:.4f}{cons})"
100
101 def gridsearch(
102 self,
103 x: ArrayLike,
104 y: ArrayLike,
105 lam: Optional[ArrayLike] = None,
106 n_splines: Optional[ArrayLike] = None,
107 gamma: float = 1.4,
108 n_jobs: int = 1,
109 weights: Optional[ArrayLike] = None,
110 ) -> "LinGAM":
111 """Select optimal (lam, n_splines) by minimising GCV."""
112 nsp = [np.atleast_1d(n_splines)] if n_splines is not None else None
113 lam_g = np.atleast_1d(lam) if lam is not None else None
114 super().gridsearch(
115 x,
116 y,
117 lam=lam_g,
118 n_splines=nsp,
119 gamma=gamma,
120 n_jobs=n_jobs,
121 weights=weights,
122 )
123 return self
124
125 def plot_decomposition(
126 self,
127 x: ArrayLike,
128 y: ArrayLike,
129 n_grid: int = 300,
130 figsize: Tuple[int, int] = (16, 5),
131 ) -> Tuple:
132 cons_str = f", constraint='{self.constraint}'" if self.constraint else ""
133 return self._plot_lingam_decomposition(
134 x,
135 y,
136 k_smooth=self.n_coefs,
137 prefix="B",
138 basis_title=f"B-Spline Basis (K_eff={self.n_coefs})",
139 fit_title="Weighted Basis Functions + Fit",
140 suptitle=f"GAM Decomposition (K={self.n_splines}{cons_str}, lam={self.lam:.4f})",
141 knots=self.knots_,
142 show_weighted=True,
143 n_grid=n_grid,
144 figsize=figsize,
145 )
magli.eilers._gam¶
Multi-term GAM with formula interface, grid search, and diagnostics.
- class magli.eilers._gam.GAMCore(formula: str, spline_order: int = 3, fit_intercept: bool = True)[source]¶
Multi-term Generalized Additive Model with formula interface.
- gridsearch(x: ArrayLike, y: ArrayLike, lam: List[ArrayLike] | ArrayLike | None = None, n_splines: List[ArrayLike] | ArrayLike | None = None, gamma: float = 1.4, search_terms: List[str | int] | str | None = None, n_jobs: int = 1, weights: ArrayLike | None = None, **kwargs) GAMCore[source]¶
Select optimal hyper-parameters via GCV grid search.
1"""Multi-term GAM with formula interface, grid search, and diagnostics."""
2
3from __future__ import annotations
4
5import os
6from concurrent.futures import ThreadPoolExecutor
7from itertools import product as cartesian_product
8from typing import Any, Dict, List, Optional, Set, Tuple, Union
9
10import numpy as np
11from numpy.typing import ArrayLike
12from scipy.linalg import block_diag, cho_factor, cho_solve
13
14from .._base import _GAMBase, _validate_weights
15from .._formula import parse_formula
16from ._solver import solve_pirls
17from ._terms import (
18 _FactorTerm,
19 _LinearTerm,
20 _SplineTerm,
21 _TensorTerm,
22 _Term,
23 instantiate_terms,
24)
25
26
27class GAMCore(_GAMBase):
28 """Multi-term Generalized Additive Model with formula interface."""
29
30 _scale_is_var = False
31
32 def __init__(
33 self,
34 formula: str,
35 spline_order: int = 3,
36 fit_intercept: bool = True,
37 ) -> None:
38 self.formula = formula
39 self.spline_order = spline_order
40 self.fit_intercept = fit_intercept
41
42 configs = parse_formula(formula)
43 self._terms: List[_Term] = instantiate_terms(configs, spline_order)
44 self._term_configs = configs
45
46 self._coef_slices: List[Tuple[int, int]] = []
47 self.coef_: Optional[np.ndarray] = None
48 self.statistics_: Optional[Dict[str, Any]] = None
49 self._design_matrix_: Optional[np.ndarray] = None
50
51 def __repr__(self) -> str:
52 status = "fitted" if self.is_fitted else "unfitted"
53 return f"GAMCore(formula='{self.formula}') [{status}]"
54
55 def _build_penalty(self) -> np.ndarray:
56 blocks = [t.build_penalty() for t in self._terms]
57 if self.fit_intercept:
58 blocks.append(np.zeros((1, 1)))
59 return blocks[0] if len(blocks) == 1 else block_diag(*blocks)
60
61 def fit(
62 self,
63 x: ArrayLike,
64 y: ArrayLike,
65 weights: Optional[ArrayLike] = None,
66 ) -> "GAMCore":
67 """Fit the model."""
68 x_arr = np.asarray(x)
69 y_arr = np.asarray(y).ravel()
70 if not np.all(np.isfinite(x_arr)) or not np.all(np.isfinite(y_arr)):
71 raise ValueError("Input data contains NaN or infinity.")
72 if x_arr.ndim == 1:
73 x_arr = x_arr.reshape(-1, 1)
74 w_arr = _validate_weights(weights, len(y_arr))
75 self._compile(x_arr)
76 self._fit_internal(x_arr, y_arr, weights=w_arr)
77 return self
78
79 def _fit_internal(
80 self,
81 x: np.ndarray,
82 y: np.ndarray,
83 weights: Optional[np.ndarray] = None,
84 gamma: float = 1.4,
85 ) -> None:
86 n = len(x)
87 X = self._build_model_matrix(x)
88 self._design_matrix_ = X
89 P = self._build_penalty()
90 self._compute_coef_slices()
91
92 if weights is not None:
93 w_sqrt = np.sqrt(weights)
94 X_w, y_w = X * w_sqrt[:, None], y * w_sqrt
95 null_deviance = float(
96 np.sum(weights * (y - np.average(y, weights=weights)) ** 2)
97 )
98 log_w_sum = float(np.sum(np.log(weights[weights > 0])))
99 else:
100 X_w, y_w = X, y
101 null_deviance = float(np.sum((y - np.mean(y)) ** 2))
102 log_w_sum = 0.0
103
104 if any(
105 isinstance(t, _SplineTerm) and t.has_shape_constraint
106 for t in self._terms
107 ):
108 coef, U1, B_solve = _solve_constrained(
109 X_w, y_w, P, self._terms, self._coef_slices, self.fit_intercept
110 )
111 else:
112 coef, U1, B_solve = solve_pirls(X_w, y_w, P)
113
114 rss = float(np.sum((y_w - X_w @ coef) ** 2))
115 self.coef_ = coef
116 edof_per_coef = (
117 np.sum(U1**2, axis=1) if n >= X.shape[1] else np.sum(U1**2, axis=0)
118 )
119 edof = float(np.sum(edof_per_coef))
120
121 scale = np.sqrt(rss / max(n - edof, 1))
122 cov = (B_solve @ B_solve.T) * scale**2
123
124 self.statistics_ = {
125 "edof": edof,
126 "edof_per_coef": edof_per_coef,
127 "scale": scale,
128 "cov": cov,
129 "se": np.sqrt(np.diag(cov)),
130 "rss": rss,
131 "null_deviance": null_deviance,
132 "n_samples": n,
133 "GCV": (n * rss) / max(n - gamma * edof, 1) ** 2,
134 "log_weights_sum": log_w_sum,
135 }
136
137 def gridsearch(
138 self,
139 x: ArrayLike,
140 y: ArrayLike,
141 lam: Optional[Union[List[ArrayLike], ArrayLike]] = None,
142 n_splines: Optional[Union[List[ArrayLike], ArrayLike]] = None,
143 gamma: float = 1.4,
144 search_terms: Optional[Union[List[Union[int, str]], str]] = None,
145 n_jobs: int = 1,
146 weights: Optional[ArrayLike] = None,
147 **kwargs,
148 ) -> "GAMCore":
149 """Select optimal hyper-parameters via GCV grid search."""
150 if lam is None and "lam_grids" in kwargs:
151 lam = kwargs.pop("lam_grids")
152 if n_splines is None and "n_splines_grids" in kwargs:
153 n_splines = kwargs.pop("n_splines_grids")
154 if kwargs:
155 raise ValueError(
156 f"Unrecognised kwargs for gridsearch(): {list(kwargs.keys())}"
157 )
158
159 x_arr = np.asarray(x)
160 y_arr = np.asarray(y).ravel()
161 if x_arr.ndim == 1:
162 x_arr = x_arr.reshape(-1, 1)
163 w_arr = _validate_weights(weights, len(y_arr))
164
165 return self._gridsearch(
166 x_arr,
167 y_arr,
168 lam=lam,
169 n_splines=n_splines,
170 gamma=gamma,
171 search_terms=search_terms,
172 n_jobs=n_jobs,
173 weights=w_arr,
174 )
175
176 def _resolve_search_terms(
177 self,
178 search_terms: Optional[Union[List[Union[int, str]], str]],
179 ) -> Set[int]:
180 if search_terms is None:
181 penalised = {
182 i
183 for i, t in enumerate(self._terms)
184 if not isinstance(t, _LinearTerm)
185 }
186 return penalised if penalised else set(range(len(self._terms)))
187 if isinstance(search_terms, str):
188 if search_terms == "all":
189 return set(range(len(self._terms)))
190 return self._indices_for_type(search_terms)
191
192 indices: Set[int] = set()
193 for item in search_terms:
194 if isinstance(item, int):
195 if item < 0 or item >= len(self._terms):
196 raise IndexError(
197 f"Term index {item} out of range (0..{len(self._terms) - 1})."
198 )
199 indices.add(item)
200 elif isinstance(item, str):
201 indices.update(self._indices_for_type(item))
202 else:
203 raise TypeError(
204 f"Invalid search_terms element: {item!r}. Expected int or str."
205 )
206 return indices
207
208 def _indices_for_type(self, kind: str) -> Set[int]:
209 type_map = {
210 "s": _SplineTerm,
211 "te": _TensorTerm,
212 "f": _FactorTerm,
213 "l": _LinearTerm,
214 }
215 kind_lower = kind.lower()
216 if kind_lower not in type_map:
217 raise ValueError(
218 f"Unknown term type '{kind}'. Expected s, te, f, or l."
219 )
220 cls = type_map[kind_lower]
221 return {i for i, t in enumerate(self._terms) if isinstance(t, cls)}
222
223 def _expand_partial_grids(
224 self,
225 user_grids: list,
226 search_idx: Set[int],
227 grid_type: str,
228 ) -> list:
229 full = []
230 user_idx = 0
231 for i, t in enumerate(self._terms):
232 if grid_type == "lam":
233 n_slots = t.d if isinstance(t, _TensorTerm) else 1
234 else:
235 n_slots = (
236 t.d
237 if isinstance(t, _TensorTerm)
238 else 1
239 if isinstance(t, _SplineTerm)
240 else 0
241 )
242 if n_slots == 0:
243 continue
244 if i in search_idx:
245 full.extend(user_grids[user_idx : user_idx + n_slots])
246 user_idx += n_slots
247 elif grid_type == "lam":
248 vals = t.lam_per if isinstance(t, _TensorTerm) else [t.lam]
249 full.extend([np.array([v]) for v in vals])
250 else:
251 vals = (
252 t.n_splines_per
253 if isinstance(t, _TensorTerm)
254 else [t.n_splines]
255 )
256 full.extend([np.array([v]) for v in vals])
257 return full
258
259 def _build_combos(
260 self,
261 x: np.ndarray,
262 lam_grids: Optional[Union[List[ArrayLike], ArrayLike]],
263 n_splines_grids: Optional[Union[List[ArrayLike], ArrayLike]],
264 search_terms: Optional[Union[List[Union[int, str]], str]] = None,
265 ) -> Tuple[List[Any], List[Tuple[float, ...]]]:
266 search_idx = self._resolve_search_terms(search_terms)
267
268 if n_splines_grids is None:
269 k_combo = [
270 int(t.n_splines)
271 if isinstance(t, _SplineTerm)
272 else tuple(int(k) for k in t.n_splines_per)
273 if isinstance(t, _TensorTerm)
274 else None
275 for t in self._terms
276 ]
277 all_k_combos = [tuple(k_combo)]
278 else:
279 n_nsp_slots = sum(
280 t.d if isinstance(t, _TensorTerm) else 1
281 for t in self._terms
282 if isinstance(t, (_SplineTerm, _TensorTerm))
283 )
284 n_nsp_searched = sum(
285 t.d if isinstance(t, _TensorTerm) else 1
286 for i, t in enumerate(self._terms)
287 if i in search_idx and isinstance(t, (_SplineTerm, _TensorTerm))
288 )
289
290 if not isinstance(n_splines_grids, list):
291 n_splines_grids = [np.atleast_1d(n_splines_grids)] * max(
292 n_nsp_slots, 1
293 )
294 elif (
295 len(n_splines_grids) == n_nsp_searched
296 and n_nsp_searched != n_nsp_slots
297 ):
298 n_splines_grids = self._expand_partial_grids(
299 n_splines_grids, search_idx, "nsp"
300 )
301 elif len(n_splines_grids) == 1 and n_nsp_slots > 1:
302 n_splines_grids = n_splines_grids * n_nsp_slots
303 elif len(n_splines_grids) != n_nsp_slots:
304 raise ValueError(
305 f"n_splines_grids has {len(n_splines_grids)} elements, "
306 f"but model requires {n_nsp_slots} n_splines slots "
307 f"({self.formula})."
308 )
309
310 min_nsp = self.spline_order + 1
311 for i, grid in enumerate(n_splines_grids):
312 for v in np.atleast_1d(grid):
313 if int(v) <= self.spline_order:
314 raise ValueError(
315 f"n_splines_grids[{i}] has value {int(v)}, "
316 f"must be > spline_order ({self.spline_order}). "
317 f"Use at least {min_nsp}."
318 )
319
320 nsp_idx = 0
321 k_combo_elements: List[List] = []
322 for t in self._terms:
323 if isinstance(t, _SplineTerm):
324 k_combo_elements.append(
325 [int(v) for v in np.atleast_1d(n_splines_grids[nsp_idx])]
326 )
327 nsp_idx += 1
328 elif isinstance(t, _TensorTerm):
329 grids = [
330 np.atleast_1d(g)
331 for g in n_splines_grids[nsp_idx : nsp_idx + t.d]
332 ]
333 k_combo_elements.append(
334 [tuple(int(v) for v in c) for c in cartesian_product(*grids)]
335 )
336 nsp_idx += t.d
337 else:
338 k_combo_elements.append([None])
339
340 all_k_combos = list(cartesian_product(*k_combo_elements))
341
342 slot_term_map: List[Tuple[int, int]] = []
343 default_fixed: List[float] = []
344 for ti, t in enumerate(self._terms):
345 if isinstance(t, _TensorTerm):
346 for d_i in range(t.d):
347 slot_term_map.append((ti, d_i))
348 default_fixed.append(float(t.lam_per[d_i]))
349 else:
350 slot_term_map.append((ti, 0))
351 default_fixed.append(float(t.lam))
352
353 total_slots = len(slot_term_map)
354 searched_slot_indices = [
355 si for si, (ti, _) in enumerate(slot_term_map) if ti in search_idx
356 ]
357
358 if lam_grids is None:
359 lam_grids = np.logspace(-3, 3, 11)
360
361 if isinstance(lam_grids, np.ndarray) and lam_grids.ndim == 2:
362 n_cols = lam_grids.shape[1]
363 if n_cols == len(searched_slot_indices):
364 lam_combos = []
365 for row in lam_grids:
366 combo = list(default_fixed)
367 for si, val in zip(searched_slot_indices, row):
368 combo[si] = float(val)
369 lam_combos.append(tuple(combo))
370 elif n_cols == total_slots:
371 lam_combos = [tuple(float(v) for v in row) for row in lam_grids]
372 else:
373 raise ValueError(
374 f"2D lam_grids has {n_cols} columns, expected {len(searched_slot_indices)} "
375 f"for searched slots or {total_slots} for all slots."
376 )
377 elif isinstance(lam_grids, list) and len(lam_grids) > 1:
378 n_searched = len(searched_slot_indices)
379 if len(lam_grids) == n_searched and n_searched != total_slots:
380 lam_grids = self._expand_partial_grids(
381 lam_grids, search_idx, "lam"
382 )
383 elif len(lam_grids) != total_slots:
384 raise ValueError(
385 f"lam_grids has {len(lam_grids)} elements, "
386 f"but model requires {total_slots} lam slots "
387 f"({self.formula})."
388 )
389 raw_combos = list(
390 cartesian_product(*[np.atleast_1d(g) for g in lam_grids])
391 )
392 final_combos = []
393 for c in raw_combos:
394 combo = list(c)
395 for si, (ti, _) in enumerate(slot_term_map):
396 if ti not in search_idx:
397 combo[si] = default_fixed[si]
398 final_combos.append(tuple(float(v) for v in combo))
399 lam_combos = list(dict.fromkeys(final_combos))
400 else:
401 sub = lam_grids[0] if isinstance(lam_grids, list) else lam_grids
402 lam_combos = []
403 for val in np.atleast_1d(sub):
404 combo = list(default_fixed)
405 for si in searched_slot_indices:
406 combo[si] = float(val)
407 lam_combos.append(tuple(combo))
408
409 self._compile(x)
410 return all_k_combos, lam_combos
411
412 def _gridsearch(
413 self,
414 x: np.ndarray,
415 y: np.ndarray,
416 lam: Optional[Union[List[ArrayLike], ArrayLike]] = None,
417 n_splines: Optional[Union[List[ArrayLike], ArrayLike]] = None,
418 gamma: float = 1.4,
419 search_terms: Optional[Union[List[Union[int, str]], str]] = None,
420 n_jobs: int = 1,
421 weights: Optional[np.ndarray] = None,
422 ) -> "GAMCore":
423 all_k_combos, lam_combos = self._build_combos(
424 x, lam, n_splines, search_terms
425 )
426 n = len(y)
427 w_sqrt = np.sqrt(weights) if weights is not None else None
428 y_w = y * w_sqrt if w_sqrt is not None else y
429 yty = float(np.dot(y_w, y_w))
430
431 def _precompute_k(k_combo):
432 configs_copy = [
433 {**cfg, "kwargs": dict(cfg["kwargs"])}
434 for cfg in self._term_configs
435 ]
436 for ci, cfg in enumerate(configs_copy):
437 k_val = k_combo[ci]
438 if cfg["type"] == "s" and k_val is not None:
439 cfg["kwargs"]["n_splines"] = k_val
440 elif cfg["type"] == "te" and k_val is not None:
441 cfg["kwargs"]["n_splines"] = (
442 list(k_val) if len(k_val) > 1 else k_val[0]
443 )
444 terms_copy = instantiate_terms(configs_copy, self.spline_order)
445 for t in terms_copy:
446 t.compile(x)
447 X_k = _build_matrix_from_terms(terms_copy, x, self.fit_intercept)
448 if w_sqrt is not None:
449 X_k = X_k * w_sqrt[:, None]
450 Q, R = np.linalg.qr(X_k)
451 qty = Q.T @ y_w
452 RtR = R.T @ R
453 Rt_qty = R.T @ qty
454 ss_orth = max(0.0, yty - float(np.dot(qty, qty)))
455 bases_list, m = _precompute_penalty_bases(
456 terms_copy, self.fit_intercept
457 )
458 return R, RtR, Rt_qty, qty, ss_orth, bases_list, m, terms_copy, X_k
459
460 n_k = len(all_k_combos)
461 workers = os.cpu_count() or 4 if n_jobs == -1 else max(1, int(n_jobs))
462 if n_jobs != 1 and n_k > 1:
463 with ThreadPoolExecutor(max_workers=min(workers, n_k)) as ex:
464 precomputed = list(ex.map(_precompute_k, all_k_combos))
465 else:
466 precomputed = [_precompute_k(c) for c in all_k_combos]
467
468 def _eval_k_combo(k_idx):
469 R, RtR, Rt_qty, qty, ss_orth, bases_list, m, terms_k, X_k = (
470 precomputed[k_idx]
471 )
472 has_shape = any(
473 isinstance(t, _SplineTerm) and t.has_shape_constraint
474 for t in terms_k
475 )
476 slices_k = [b[1] for b in bases_list]
477 P_buf = np.zeros((m, m))
478 B_buf = np.empty((m, m))
479 best_gcv_k = np.inf
480 best_l_idx = 0
481
482 for l_idx, l_combo in enumerate(lam_combos):
483 _assemble_penalty_from_bases(bases_list, l_combo, m, out=P_buf)
484 if has_shape:
485 coef, U1, _ = _solve_constrained(
486 X_k, y_w, P_buf, terms_k, slices_k, self.fit_intercept
487 )
488 edf = float(np.sum(U1**2))
489 rss = float(np.sum((y_w - X_k @ coef) ** 2))
490 else:
491 np.add(RtR, P_buf, out=B_buf)
492 B_buf.flat[:: m + 1] += 1e-8
493 try:
494 c_chol, lower = cho_factor(
495 B_buf, overwrite_a=True, check_finite=False
496 )
497 coef = cho_solve(
498 (c_chol, lower), Rt_qty, check_finite=False
499 )
500 edf = float(
501 np.trace(
502 cho_solve(
503 (c_chol, lower), RtR, check_finite=False
504 )
505 )
506 )
507 except Exception:
508 B_fallback = RtR + P_buf
509 B_fallback.flat[:: m + 1] += 1e-8
510 coef, _, _, _ = np.linalg.lstsq(
511 B_fallback, Rt_qty, rcond=1e-10
512 )
513 edf = float(
514 np.trace(
515 np.linalg.pinv(B_fallback, rcond=1e-10) @ RtR
516 )
517 )
518 r_diff = qty - R @ coef
519 rss = max(0.0, ss_orth + float(np.dot(r_diff, r_diff)))
520
521 gcv = (n * rss) / max(n - gamma * edf, 1) ** 2
522 if gcv < best_gcv_k:
523 best_gcv_k = gcv
524 best_l_idx = l_idx
525
526 return best_gcv_k, k_idx, best_l_idx
527
528 if n_jobs != 1 and n_k > 1:
529 with ThreadPoolExecutor(max_workers=min(workers, n_k)) as ex:
530 k_results = list(ex.map(_eval_k_combo, range(n_k)))
531 else:
532 k_results = [_eval_k_combo(i) for i in range(n_k)]
533
534 _, best_k_idx, best_l_idx = min(k_results, key=lambda r: r[0])
535 self._terms = precomputed[best_k_idx][-2]
536 _set_lams_from_combo(self._terms, lam_combos[best_l_idx])
537 self._compile(x)
538 self._fit_internal(x, y, weights=weights, gamma=gamma)
539 return self
540
541 def plot_decomposition(
542 self,
543 x: ArrayLike,
544 y: ArrayLike,
545 n_grid: int = 200,
546 extrap_frac: float = 0.15,
547 figsize: Optional[Tuple[int, int]] = None,
548 feature_names: Optional[List[str]] = None,
549 ) -> Tuple:
550 return self._plot_gamcore_decomposition(
551 x,
552 y,
553 n_grid=n_grid,
554 extrap_frac=extrap_frac,
555 figsize=figsize,
556 feature_names=feature_names,
557 )
558
559
560def _solve_constrained(
561 X_w: np.ndarray,
562 y_w: np.ndarray,
563 P: np.ndarray,
564 terms: List[_Term],
565 coef_slices: List[Tuple[int, int]],
566 fit_intercept: bool,
567 max_iter: int = 100,
568 tol: float = 1e-4,
569) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
570 m = X_w.shape[1]
571 coef = np.ones(m) * np.sqrt(np.finfo(float).eps)
572 for _ in range(max_iter):
573 blocks = [
574 t.build_constraint_penalty(coef[s:e])
575 if isinstance(t, _SplineTerm)
576 else np.zeros((t.n_coefs, t.n_coefs))
577 for t, (s, e) in zip(terms, coef_slices)
578 ]
579 if fit_intercept:
580 blocks.append(np.zeros((1, 1)))
581 C = blocks[0] if len(blocks) == 1 else block_diag(*blocks)
582 coef_new, U1, B_solve = solve_pirls(X_w, y_w, P + C)
583 if np.linalg.norm(coef - coef_new) < tol:
584 return coef_new, U1, B_solve
585 coef = coef_new
586 return coef, U1, B_solve
587
588
589def _set_lams_from_combo(terms: List[_Term], l_combo: tuple) -> None:
590 lam_idx = 0
591 for t in terms:
592 if isinstance(t, _TensorTerm):
593 t.lam_per = list(l_combo[lam_idx : lam_idx + t.d])
594 lam_idx += t.d
595 else:
596 t.lam = float(l_combo[lam_idx])
597 lam_idx += 1
598
599
600def _assemble_penalty_from_bases(
601 bases_list: list,
602 l_combo: tuple,
603 m: int,
604 out: Optional[np.ndarray] = None,
605) -> np.ndarray:
606 P = np.zeros((m, m)) if out is None else out
607 if out is not None:
608 P[:] = 0
609 lam_idx = 0
610 for base_mats, (start, end), n_slots in bases_list:
611 for k in range(n_slots):
612 P[start:end, start:end] += float(l_combo[lam_idx]) * base_mats[k]
613 lam_idx += 1
614 return P
615
616
617def _precompute_penalty_bases(
618 terms: List[_Term],
619 fit_intercept: bool,
620) -> Tuple[list, int]:
621 bases_list = []
622 off = 0
623 for t in terms:
624 pb = t.penalty_bases()
625 bases_list.append(
626 ([b[0] for b in pb], (off, off + t.n_coefs), sum(b[1] for b in pb))
627 )
628 off += t.n_coefs
629 return bases_list, off + (1 if fit_intercept else 0)
630
631
632def _build_matrix_from_terms(
633 terms: List[_Term],
634 x: np.ndarray,
635 fit_intercept: bool,
636) -> np.ndarray:
637 blocks = [t.build_columns(x) for t in terms]
638 if fit_intercept:
639 blocks.append(np.ones((len(x), 1)))
640 return np.hstack(blocks)
Wood Engine (magli.wood)¶
Python implementation of Simon Wood’s GAM methods (Wood 2011, 2017) with multiple basis types, null-space constraint absorption, and continuous REML / GCV smoothing optimization.
Classes¶
Modules¶
magli.wood._smooths¶
Basis constructors and constraint absorption following Wood (2011, 2017).
Supported basis types: - ‘tp’: Thin plate regression splines (default in Wood backend) - ‘ts’: Thin plate regression splines with shrinkage - ‘cr’: Cubic regression splines (knots at quantiles, natural boundaries) - ‘cs’: Cubic regression splines with shrinkage - ‘cc’: Cyclic cubic regression splines (Lancaster & Salkauskas) - ‘ps’: P-splines (B-splines with discrete difference penalties) - ‘cp’: Cyclic P-splines - ‘bs’: B-splines with integrated derivative penalties - ‘re’: Random effects / ridge penalties
- magli.wood._smooths.build_cr_basis(x: ndarray, k: int = 10, knots: ndarray | None = None, shrink: bool = False) Tuple[ndarray, List[ndarray], Dict][source]¶
Cubic regression spline basis (‘cr’) or with shrinkage (‘cs’).
- magli.wood._smooths.build_cc_basis(x: ndarray, k: int = 10, knots: ndarray | None = None) Tuple[ndarray, List[ndarray], Dict][source]¶
Cyclic cubic regression spline (‘cc’) based on Lancaster & Salkauskas (1986).
- magli.wood._smooths.build_ps_basis(x: ndarray, k: int = 10, knots: ndarray | None = None, m: Tuple[int, int] = (2, 2)) Tuple[ndarray, List[ndarray], Dict][source]¶
P-spline basis (‘ps’) following Eilers & Marx and Wood (2017).
- magli.wood._smooths.build_cp_basis(x: ndarray, k: int = 10, knots: ndarray | None = None, m: Tuple[int, int] = (2, 2)) Tuple[ndarray, List[ndarray], Dict][source]¶
Cyclic P-spline basis (‘cp’).
- magli.wood._smooths.build_bs_basis(x: ndarray, k: int = 10, knots: ndarray | None = None, m: Tuple[int, int] = (3, 2)) Tuple[ndarray, List[ndarray], Dict][source]¶
B-spline basis (‘bs’) with integrated squared derivative penalty.
- magli.wood._smooths.build_tp_basis(x: ndarray, k: int = 10, m: int = 2, shrink: bool = False) Tuple[ndarray, List[ndarray], Dict][source]¶
Thin plate regression spline basis (‘tp’ / ‘ts’) following Wood (2003).
- magli.wood._smooths.predict_smooth_basis(x: ndarray, meta: Dict[str, Any]) ndarray[source]¶
Evaluate basis functions on new covariate values x using stored metadata.
- magli.wood._smooths.build_smooth_basis(x: ndarray, bs: str = 'tp', k: int = 10, knots: ndarray | None = None, m: int | Tuple[int, int] | None = None) Tuple[ndarray, List[ndarray], Dict][source]¶
Dispatch to the requested basis constructor.
- magli.wood._smooths.rescale_penalty(X: ndarray, S_list: List[ndarray]) Tuple[List[ndarray], List[float]][source]¶
Rescale penalty matrix by norm(X, ‘I’)^2.
- magli.wood._smooths.absorb_sum_to_zero_constraint(X: ndarray, S_list: List[ndarray]) Tuple[ndarray, List[ndarray], ndarray][source]¶
Absorb sample sum-to-zero constraint C beta = 0 where C = 1^T X.
1"""Basis constructors and constraint absorption following Wood (2011, 2017).
2
3Supported basis types:
4- 'tp': Thin plate regression splines (default in Wood backend)
5- 'ts': Thin plate regression splines with shrinkage
6- 'cr': Cubic regression splines (knots at quantiles, natural boundaries)
7- 'cs': Cubic regression splines with shrinkage
8- 'cc': Cyclic cubic regression splines (Lancaster & Salkauskas)
9- 'ps': P-splines (B-splines with discrete difference penalties)
10- 'cp': Cyclic P-splines
11- 'bs': B-splines with integrated derivative penalties
12- 're': Random effects / ridge penalties
13"""
14
15from __future__ import annotations
16
17from typing import Any, Dict, List, Optional, Tuple, Union
18
19import numpy as np
20from scipy import linalg
21from scipy.interpolate import BSpline
22
23
24def _place_knots(x: np.ndarray, nk: int) -> np.ndarray:
25 """R's place.knots: evenly spaced knot placement throughout unique x."""
26 xu = np.sort(np.unique(x))
27 n = len(xu)
28 if nk > n:
29 raise ValueError("More knots than unique data values is not allowed.")
30 if nk < 2:
31 raise ValueError("Too few knots.")
32 if nk == 2:
33 return np.array([xu[0], xu[-1]], dtype=float)
34
35 delta = (n - 1) / (nk - 1)
36 indices = np.arange(1, nk - 1)
37 lbi = np.floor(delta * indices).astype(int)
38 frac = delta * indices - lbi
39
40 knots = np.zeros(nk, dtype=float)
41 knots[0] = xu[0]
42 knots[-1] = xu[-1]
43 knots[1:-1] = xu[lbi] * (1.0 - frac) + xu[1:][lbi] * frac
44 return knots
45
46
47def _get_cr_FS(xk: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
48 """Compute wiggliness penalty S and matrix F mapping knot values to second derivs."""
49 n = len(xk)
50 h = np.diff(xk)
51 n2 = n - 2
52
53 D = np.zeros((n2, n), dtype=float)
54 idx = np.arange(n2)
55 D[idx, idx] = 1.0 / h[:-1]
56 D[idx, idx + 1] = -1.0 / h[:-1] - 1.0 / h[1:]
57 D[idx, idx + 2] = 1.0 / h[1:]
58
59 ab = np.zeros((2, n2), dtype=float)
60 ab[0, :] = (h[:-1] + h[1:]) / 3.0
61 ab[1, : n2 - 1] = h[1:-1] / 6.0
62 B_inv_D = linalg.solveh_banded(ab, D, lower=True)
63
64 F = np.zeros((n, n), dtype=float)
65 F[1:-1, :] = B_inv_D
66
67 S = D.T @ B_inv_D
68 return F, 0.5 * (S + S.T)
69
70
71def _eval_cr_basis(x: np.ndarray, xk: np.ndarray, F: np.ndarray) -> np.ndarray:
72 """Evaluate cubic regression spline basis at points x with given knots and F."""
73 n = len(x)
74 nk = len(xk)
75 X = np.zeros((n, nk), dtype=float)
76 kmin, kmax = xk[0], xk[-1]
77
78 m_low = x < kmin
79 m_high = x > kmax
80 m_mid = ~(m_low | m_high)
81
82 if np.any(m_low):
83 h = xk[1] - kmin
84 xik = x[m_low] - kmin
85 X[m_low, :] = (-xik * (h / 3.0))[:, None] * F[0, :] + (
86 -xik * (h / 6.0)
87 )[:, None] * F[1, :]
88 X[m_low, 0] += 1.0 - xik / h
89 X[m_low, 1] += xik / h
90
91 if np.any(m_high):
92 j = nk - 1
93 h = kmax - xk[j - 1]
94 xik = x[m_high] - kmax
95 X[m_high, :] = (xik * (h / 6.0))[:, None] * F[j - 1, :] + (
96 xik * (h / 3.0)
97 )[:, None] * F[j, :]
98 X[m_high, nk - 2] += -xik / h
99 X[m_high, nk - 1] += 1.0 + xik / h
100
101 if np.any(m_mid):
102 x_mid = x[m_mid]
103 j = np.clip(np.searchsorted(xk, x_mid, side="right") - 1, 0, nk - 2)
104 xj, xj1 = xk[j], xk[j + 1]
105 h = xj1 - xj
106 ajm, ajp = xj1 - x_mid, x_mid - xj
107 cjm = (ajm * (ajm * ajm / h - h)) / 6.0
108 cjp = (ajp * (ajp * ajp / h - h)) / 6.0
109
110 idx_mid = np.where(m_mid)[0]
111 X[idx_mid, :] = cjm[:, None] * F[j, :] + cjp[:, None] * F[j + 1, :]
112 np.add.at(X, (idx_mid, j), ajm / h)
113 np.add.at(X, (idx_mid, j + 1), ajp / h)
114
115 return X
116
117
118def build_cr_basis(
119 x: np.ndarray,
120 k: int = 10,
121 knots: Optional[np.ndarray] = None,
122 shrink: bool = False,
123) -> Tuple[np.ndarray, List[np.ndarray], Dict]:
124 """Cubic regression spline basis ('cr') or with shrinkage ('cs')."""
125 x = np.asarray(x, dtype=float).ravel()
126 k = max(k, 3)
127 if knots is None:
128 xu = np.sort(np.unique(x))
129 if len(xu) < k:
130 raise ValueError(
131 f"Insufficient unique covariate values ({len(xu)}) for k={k} knots."
132 )
133 xk = np.quantile(xu, np.linspace(0.0, 1.0, k))
134 else:
135 xk = np.asarray(knots, dtype=float).ravel()
136 if len(xk) != k:
137 raise ValueError(
138 f"Number of supplied knots ({len(xk)}) != k ({k}) for cr smooth."
139 )
140
141 F, S = _get_cr_FS(xk)
142 X = _eval_cr_basis(x, xk, F)
143
144 if shrink:
145 es_val, es_vec = np.linalg.eigh(S)
146 es_val[1] = es_val[2] * 0.1
147 es_val[0] = es_val[1] * 0.1
148 S = es_vec @ (es_val[:, None] * es_vec.T)
149 S = 0.5 * (S + S.T)
150
151 meta = {
152 "bs": "cs" if shrink else "cr",
153 "xp": xk,
154 "F": F,
155 "rank": k if shrink else k - 2,
156 "null_space_dim": 0 if shrink else 2,
157 }
158 return X, [S], meta
159
160
161def _cwrap(x0: float, x1: float, x: np.ndarray) -> np.ndarray:
162 """Map x onto [x0, x1] cyclically."""
163 h = x1 - x0
164 x_out = x.copy()
165 high = x_out > x1
166 if np.any(high):
167 x_out[high] = x0 + (x_out[high] - x1) % h
168 low = x_out < x0
169 if np.any(low):
170 x_out[low] = x1 - (x0 - x_out[low]) % h
171 return x_out
172
173
174def build_cc_basis(
175 x: np.ndarray,
176 k: int = 10,
177 knots: Optional[np.ndarray] = None,
178) -> Tuple[np.ndarray, List[np.ndarray], Dict]:
179 """Cyclic cubic regression spline ('cc') based on Lancaster & Salkauskas (1986)."""
180 x = np.asarray(x, dtype=float).ravel()
181 k = max(k, 4)
182 if knots is None:
183 xk = _place_knots(x, k)
184 elif len(knots) == 2:
185 xk = _place_knots(np.concatenate([knots, x]), k)
186 else:
187 xk = np.asarray(knots, dtype=float).ravel()
188 if len(xk) != k:
189 raise ValueError(
190 f"Number of supplied knots ({len(xk)}) != k ({k}) for cc smooth."
191 )
192
193 h = np.diff(xk)
194 n = len(xk) - 1
195 idx = np.arange(n)
196 prev_i = (idx - 1) % n
197 next_i = (idx + 1) % n
198
199 B = np.zeros((n, n), dtype=float)
200 D = np.zeros((n, n), dtype=float)
201 B[idx, prev_i] = h[prev_i] / 6.0
202 B[idx, idx] = (h[prev_i] + h[idx]) / 3.0
203 B[idx, next_i] = h[idx] / 6.0
204 D[idx, prev_i] = 1.0 / h[prev_i]
205 D[idx, idx] = -(1.0 / h[prev_i] + 1.0 / h[idx])
206 D[idx, next_i] = 1.0 / h[idx]
207
208 BD = np.linalg.solve(B, D)
209 S = D.T @ BD
210 S = 0.5 * (S + S.T)
211 X = _eval_cc_basis(x, xk, BD)
212
213 return X, [S], {"bs": "cc", "xp": xk, "BD": BD, "rank": n - 1, "null_space_dim": 1}
214
215
216def _eval_cc_basis(x: np.ndarray, xk: np.ndarray, BD: np.ndarray) -> np.ndarray:
217 """Evaluate cyclic cubic regression spline basis at points x with given knots and BD."""
218 x_wrapped = _cwrap(xk[0], xk[-1], x)
219 nx = len(x_wrapped)
220 n = len(xk) - 1
221 h = np.diff(xk)
222
223 j = np.clip(np.searchsorted(xk, x_wrapped, side="right") - 1, 0, n - 1)
224 hj = h[j]
225 c1 = (xk[j + 1] - x_wrapped) ** 3 / (6.0 * hj)
226 c2 = (x_wrapped - xk[j]) ** 3 / (6.0 * hj)
227 c3 = (hj * (xk[j + 1] - x_wrapped)) / 6.0
228 c4 = (hj * (x_wrapped - xk[j])) / 6.0
229 c5 = (xk[j + 1] - x_wrapped) / hj
230 c6 = (x_wrapped - xk[j]) / hj
231
232 j_next = np.where(j + 1 == n, 0, j + 1)
233 X = BD[j, :] * (c1 - c3)[:, None] + BD[j_next, :] * (c2 - c4)[:, None]
234 idx = np.arange(nx)
235 np.add.at(X, (idx, j), c5)
236 np.add.at(X, (idx, j_next), c6)
237 return X
238
239
240def _spline_design(
241 knots: np.ndarray, x: np.ndarray, ord: int, outer_ok: bool = False
242) -> np.ndarray:
243 """Evaluate standard B-spline design matrix."""
244 n_bases = len(knots) - ord
245 k_deg = ord - 1
246 k_pad = np.concatenate(
247 [
248 np.repeat(knots[0] - 1.0, k_deg),
249 knots,
250 np.repeat(knots[-1] + 1.0, k_deg),
251 ]
252 )
253 c = np.zeros((len(k_pad) - ord, n_bases), dtype=float)
254 c[k_deg : k_deg + n_bases, :] = np.eye(n_bases)
255 out = BSpline(k_pad, c, k_deg, extrapolate=False)(x)
256 np.nan_to_num(out, copy=False, nan=0.0)
257 return out
258
259
260def _c_spline_des(x: np.ndarray, knots: np.ndarray, ord: int = 4) -> np.ndarray:
261 """Cyclic B-spline basis."""
262 knots = np.sort(np.asarray(knots, dtype=float).ravel())
263 nk = len(knots)
264 k1, kmax = knots[0], knots[-1]
265 x_w = _cwrap(k1, kmax, x)
266 xc = knots[nk - ord]
267 knots_ext = np.concatenate([k1 - (kmax - knots[nk - ord : nk - 1]), knots])
268 X = _spline_design(knots_ext, x_w, ord, outer_ok=True)
269 ind = x_w > xc
270 if np.any(ind):
271 X[ind] += _spline_design(
272 knots_ext, x_w[ind] - kmax + k1, ord, outer_ok=True
273 )
274 return X
275
276
277def build_ps_basis(
278 x: np.ndarray,
279 k: int = 10,
280 knots: Optional[np.ndarray] = None,
281 m: Tuple[int, int] = (2, 2),
282) -> Tuple[np.ndarray, List[np.ndarray], Dict]:
283 """P-spline basis ('ps') following Eilers & Marx and Wood (2017)."""
284 x = np.asarray(x, dtype=float).ravel()
285 m_basis, m_pen = m
286 k = max(k, m_basis + 1)
287 nk_interior = k - m_basis
288
289 if knots is None or len(knots) == 2:
290 xl, xu = (
291 (float(knots[0]), float(knots[1]))
292 if knots is not None and len(knots) == 2
293 else (float(np.min(x)), float(np.max(x)))
294 )
295 xr = xu - xl
296 xl -= xr * 0.001
297 xu += xr * 0.001
298 dx = (xu - xl) / (nk_interior - 1)
299 knots_full = np.linspace(
300 xl - dx * (m_basis + 1),
301 xu + dx * (m_basis + 1),
302 nk_interior + 2 * m_basis + 2,
303 )
304 else:
305 knots_full = np.asarray(knots, dtype=float).ravel()
306
307 order = m_basis + 2
308 X = _spline_design(knots_full, x, order, outer_ok=True)
309 D = np.diff(np.eye(k), n=m_pen, axis=0)
310 S = D.T @ D
311 S = 0.5 * (S + S.T)
312
313 return (
314 X,
315 [S],
316 {
317 "bs": "ps",
318 "knots": knots_full,
319 "ord": order,
320 "m": m,
321 "rank": k - m_pen,
322 "null_space_dim": m_pen,
323 },
324 )
325
326
327def build_cp_basis(
328 x: np.ndarray,
329 k: int = 10,
330 knots: Optional[np.ndarray] = None,
331 m: Tuple[int, int] = (2, 2),
332) -> Tuple[np.ndarray, List[np.ndarray], Dict]:
333 """Cyclic P-spline basis ('cp')."""
334 x = np.asarray(x, dtype=float).ravel()
335 m_basis, m_pen = m
336 k = max(k, m_basis, 4)
337 nk = k + 1
338
339 if knots is None or len(knots) == 2:
340 x0, x1 = (
341 (float(np.min(knots)), float(np.max(knots)))
342 if knots is not None and len(knots) == 2
343 else (float(np.min(x)), float(np.max(x)))
344 )
345 knots_seq = np.linspace(x0, x1, nk)
346 else:
347 knots_seq = np.asarray(knots, dtype=float).ravel()
348 if len(knots_seq) != nk:
349 raise ValueError(f"there should be {nk} supplied knots")
350
351 order = m_basis + 2
352 X = _c_spline_des(x, knots_seq, ord=order)
353
354 if m_pen > 0:
355 De = np.diff(np.eye(k + m_pen), n=m_pen, axis=0)
356 D = De[:, m_pen:].copy()
357 D[:, -m_pen:] += De[:, :m_pen]
358 else:
359 D = np.eye(k)
360 S = D.T @ D
361 S = 0.5 * (S + S.T)
362
363 return (
364 X,
365 [S],
366 {
367 "bs": "cp",
368 "knots": knots_seq,
369 "ord": order,
370 "m": m,
371 "rank": k - 1,
372 "null_space_dim": 1,
373 },
374 )
375
376
377def build_bs_basis(
378 x: np.ndarray,
379 k: int = 10,
380 knots: Optional[np.ndarray] = None,
381 m: Tuple[int, int] = (3, 2),
382) -> Tuple[np.ndarray, List[np.ndarray], Dict]:
383 """B-spline basis ('bs') with integrated squared derivative penalty."""
384 x = np.asarray(x, dtype=float).ravel()
385 m1, m2 = m
386 k = max(k, m1 + 1)
387 nk = k - m1 + 1
388
389 if knots is None or len(knots) == 2:
390 xl, xu = (
391 (float(knots[0]), float(knots[1]))
392 if knots is not None and len(knots) == 2
393 else (float(np.min(x)), float(np.max(x)))
394 )
395 xr = xu - xl
396 xl -= xr * 0.001
397 xu += xr * 0.001
398 dx = (xu - xl) / (nk - 1)
399 knots_full = np.linspace(xl - dx * m1, xu + dx * m1, nk + 2 * m1)
400 else:
401 knots_full = np.asarray(knots, dtype=float).ravel()
402
403 order = m1 + 1
404 X = _spline_design(knots_full, x, order, outer_ok=True)
405
406 k0 = knots_full[m1 : m1 + nk]
407 h = np.diff(k0)
408 pord = m1 - m2
409 if pord > 0:
410 xi, wi = np.polynomial.legendre.leggauss(pord + 1)
411 pts = (0.5 * (k0[:-1, None] + k0[1:, None]) + 0.5 * h[:, None] * xi[None, :]).ravel()
412 wts = (0.5 * h[:, None] * wi[None, :]).ravel()
413 else:
414 pts = 0.5 * (k0[:-1] + k0[1:])
415 wts = h
416
417 k_pad = np.concatenate(
418 [
419 np.repeat(knots_full[0] - 1.0, m1),
420 knots_full,
421 np.repeat(knots_full[-1] + 1.0, m1),
422 ]
423 )
424 c = np.zeros((len(k_pad) - order, k), dtype=float)
425 c[m1 : m1 + k, :] = np.eye(k)
426 G = BSpline(k_pad, c, m1, extrapolate=False)(pts, nu=m2)
427 np.nan_to_num(G, copy=False, nan=0.0)
428 S = G.T @ (wts[:, None] * G)
429 S = 0.5 * (S + S.T)
430
431 return (
432 X,
433 [S],
434 {
435 "bs": "bs",
436 "knots": knots_full,
437 "ord": order,
438 "m": m,
439 "rank": k - m2,
440 "null_space_dim": m2,
441 },
442 )
443
444
445def _eta_1d(r: np.ndarray) -> np.ndarray:
446 """Thin plate spline Green's function for d=1, m=2: eta(r) = (1/12) * |r|^3."""
447 out = np.abs(r)
448 out *= r * r * (1.0 / 12.0)
449 return out
450
451
452def _wood_qt(A: np.ndarray) -> np.ndarray:
453 """Right-hand Householder factorization: A Q = [0, R]."""
454 M, k = A.shape
455 Q = np.eye(k)
456 A_work = A.copy()
457 for i in range(M):
458 p = A_work[i, : k - i].copy()
459 m = np.max(np.abs(p))
460 if m > 0.0:
461 p /= m
462 lsq = np.sqrt(np.sum(p * p))
463 if p[k - i - 1] < 0.0:
464 lsq = -lsq
465 p[k - i - 1] += lsq
466 g = 1.0 / (lsq * p[k - i - 1]) if lsq != 0.0 else 0.0
467 lsq *= m
468 A_work[i + 1 :, : k - i] -= (A_work[i + 1 :, : k - i] @ p)[:, None] * (
469 g * p
470 )
471 Q[:, : k - i] -= (Q[:, : k - i] @ p)[:, None] * (g * p)
472 A_work[i, k - i - 1] = -lsq
473 A_work[i, : k - i - 1] = 0.0
474 return Q
475
476
477def build_tp_basis(
478 x: np.ndarray,
479 k: int = 10,
480 m: int = 2,
481 shrink: bool = False,
482) -> Tuple[np.ndarray, List[np.ndarray], Dict]:
483 """Thin plate regression spline basis ('tp' / 'ts') following Wood (2003)."""
484 x = np.asarray(x, dtype=float).ravel()
485 M = 2
486 k = max(k, M + 1)
487
488 x_mean = float(np.mean(x))
489 x_c = x - x_mean
490
491 xu, inv_idx = np.unique(x_c, return_inverse=True)
492 nu = len(xu)
493 if nu < k:
494 raise ValueError(
495 f"Covariate has only {nu} unique values, cannot support basis dimension k={k}."
496 )
497
498 E = _eta_1d(xu[:, None] - xu[None, :])
499 T = np.column_stack([np.ones(nu), xu])
500
501 eigvals, eigvecs = np.linalg.eigh(E)
502 top_k = np.argsort(np.abs(eigvals))[::-1][:k]
503 sort_idx = top_k[np.argsort(eigvals[top_k])[::-1]]
504 v_k, U_k = eigvals[sort_idx], eigvecs[:, sort_idx]
505
506 Q = _wood_qt(T.T @ U_k)
507 X_u = (U_k * v_k[None, :]) @ Q
508 X_u[:, k - M :] = T
509 UkZ = U_k @ Q[:, : k - M]
510
511 S_full = Q.T @ (v_k[:, None] * Q)
512 S_full[k - M :, :] = 0.0
513 S_full[:, k - M :] = 0.0
514 S = 0.5 * (S_full + S_full.T)
515
516 X = X_u[inv_idx, :]
517 col_rms = np.sqrt(np.mean(X**2, axis=0))
518 col_rms[col_rms == 0.0] = 1.0
519 X = X / col_rms[None, :]
520 S = S / (col_rms[:, None] * col_rms[None, :])
521 S = 0.5 * (S + S.T)
522
523 if shrink:
524 es_val, es_vec = np.linalg.eigh(S)
525 es_val[:M] = es_val[M] * 0.1
526 S = es_vec @ (es_val[:, None] * es_vec.T)
527 S = 0.5 * (S + S.T)
528
529 return (
530 X,
531 [S],
532 {
533 "bs": "ts" if shrink else "tp",
534 "shift": x_mean,
535 "col_rms": col_rms,
536 "Xu": xu,
537 "UkZ": UkZ,
538 "k": k,
539 "M": M,
540 "rank": k if shrink else k - M,
541 "null_space_dim": 0 if shrink else M,
542 },
543 )
544
545
546def predict_smooth_basis(x: np.ndarray, meta: Dict[str, Any]) -> np.ndarray:
547 """Evaluate basis functions on new covariate values x using stored metadata."""
548 x = np.asarray(x, dtype=float).ravel()
549 bs = meta["bs"]
550 if bs in ("tp", "ts"):
551 shift = meta["shift"]
552 X_smooth = _eta_1d((x - shift)[:, None] - meta["Xu"][None, :]) @ meta["UkZ"]
553 T_new = np.column_stack([np.ones(len(x)), x - shift])
554 return np.column_stack([X_smooth, T_new]) / meta["col_rms"][None, :]
555 if bs in ("cr", "cs"):
556 return _eval_cr_basis(x, meta["xp"], meta["F"])
557 if bs == "cc":
558 return _eval_cc_basis(x, meta["xp"], meta["BD"])
559 if bs == "cp":
560 return _c_spline_des(x, meta["knots"], meta["ord"])
561 if bs in ("ps", "bs", "b"):
562 return _spline_design(meta["knots"], x, meta["ord"], outer_ok=True)
563 if bs == "re":
564 return (x[:, None] == meta["levels"][None, :]).astype(float)
565 raise ValueError(f"Unknown basis type bs='{bs}' in predict_smooth_basis.")
566
567
568def build_smooth_basis(
569 x: np.ndarray,
570 bs: str = "tp",
571 k: int = 10,
572 knots: Optional[np.ndarray] = None,
573 m: Optional[Union[int, Tuple[int, int]]] = None,
574) -> Tuple[np.ndarray, List[np.ndarray], Dict]:
575 """Dispatch to the requested basis constructor."""
576 bs = bs.lower()
577 if bs in ("tp", "ts"):
578 return build_tp_basis(
579 x, k=k, m=m if isinstance(m, int) else 2, shrink=(bs == "ts")
580 )
581 if bs in ("cr", "cs"):
582 return build_cr_basis(x, k=k, knots=knots, shrink=(bs == "cs"))
583 if bs == "cc":
584 return build_cc_basis(x, k=k, knots=knots)
585 if bs == "ps":
586 return build_ps_basis(
587 x, k=k, knots=knots, m=m if isinstance(m, tuple) else (2, 2)
588 )
589 if bs == "cp":
590 return build_cp_basis(
591 x, k=k, knots=knots, m=m if isinstance(m, tuple) else (2, 2)
592 )
593 if bs in ("bs", "b"):
594 return build_bs_basis(
595 x, k=k, knots=knots, m=m if isinstance(m, tuple) else (3, 2)
596 )
597 if bs == "re":
598 xu, inv_idx = np.unique(x, return_inverse=True)
599 K = len(xu)
600 return (
601 np.eye(K)[inv_idx],
602 [np.eye(K)],
603 {"bs": "re", "levels": xu, "rank": K, "null_space_dim": 0},
604 )
605 raise ValueError(
606 f"Unknown basis type bs='{bs}'. Supported types: 'tp', 'ts', 'cr', 'cs', 'cc', 'ps', 'cp', 'bs', 're'."
607 )
608
609
610def rescale_penalty(
611 X: np.ndarray, S_list: List[np.ndarray]
612) -> Tuple[List[np.ndarray], List[float]]:
613 """Rescale penalty matrix by norm(X, 'I')^2."""
614 if not S_list:
615 return [], []
616 maXX = float(np.max(np.sum(np.abs(X), axis=1)) ** 2) or 1.0
617 S_scaled, scales = [], []
618 for S in S_list:
619 maS = float(np.max(np.sum(np.abs(S), axis=0)) / maXX) or 1.0
620 S_scaled.append(S / maS)
621 scales.append(maS)
622 return S_scaled, scales
623
624
625def absorb_sum_to_zero_constraint(
626 X: np.ndarray,
627 S_list: List[np.ndarray],
628) -> Tuple[np.ndarray, List[np.ndarray], np.ndarray]:
629 """Absorb sample sum-to-zero constraint C beta = 0 where C = 1^T X."""
630 C = np.sum(X, axis=0, keepdims=True)
631 Q, _ = np.linalg.qr(C.T, mode="complete")
632 Z = Q[:, 1:]
633 X_tilde = X @ Z
634 S_tilde = []
635 for S in S_list:
636 ZS = Z.T @ S @ Z
637 S_tilde.append(0.5 * (ZS + ZS.T))
638 return X_tilde, S_tilde, Z
magli.wood._optimizer¶
Penalized least squares solvers and smoothing parameter optimization for Wood (2011) GAM backend.
- magli.wood._optimizer.fit_wood_engine(X: ndarray, y: ndarray, S_list: List[ndarray], sp_fixed: List[float | None] | None = None, method: str = 'REML', gamma: float = 1.0, coef_slices: List[Tuple[int, int]] | None = None, M_p: int = 1, weights: ndarray | None = None) Dict[str, Any][source]¶
Fit GAM using Wood (2011) REML or GCV criteria and compute full statistics.
1"""Penalized least squares solvers and smoothing parameter optimization for Wood (2011) GAM backend."""
2
3from __future__ import annotations
4
5from typing import Any, Dict, List, Optional, Tuple
6
7import numpy as np
8from scipy import linalg, optimize, stats
9
10
11def fit_wood_engine(
12 X: np.ndarray,
13 y: np.ndarray,
14 S_list: List[np.ndarray],
15 sp_fixed: Optional[List[Optional[float]]] = None,
16 method: str = "REML",
17 gamma: float = 1.0,
18 coef_slices: Optional[List[Tuple[int, int]]] = None,
19 M_p: int = 1,
20 weights: Optional[np.ndarray] = None,
21) -> Dict[str, Any]:
22 """Fit GAM using Wood (2011) REML or GCV criteria and compute full statistics."""
23 n, p = X.shape
24 method = method.upper()
25 if method not in ("REML", "GCV.CP", "GCV"):
26 raise ValueError(f"Unknown method '{method}'. Choose 'REML' or 'GCV.Cp'.")
27
28 if weights is not None:
29 null_deviance = float(
30 np.sum(weights * (y - float(np.average(y, weights=weights))) ** 2)
31 )
32 log_w_sum = float(np.sum(np.log(weights[weights > 0])))
33 w_sqrt = np.sqrt(weights)
34 X = X * w_sqrt[:, None]
35 y = y * w_sqrt
36 else:
37 null_deviance = float(np.sum((y - float(np.mean(y))) ** 2))
38 log_w_sum = 0.0
39
40 if n > p:
41 R_aug = linalg.qr(
42 np.column_stack([X, y]), mode="r", check_finite=False
43 )[0]
44 R = R_aug[:p, :p]
45 qty = R_aug[:p, p]
46 ss_orth = float(R_aug[p, p] ** 2)
47 else:
48 Q, R = np.linalg.qr(X)
49 qty = Q.T @ y
50 ss_orth = max(0.0, float(y @ y) - float(qty @ qty))
51 RtR = R.T @ R
52 Xty = R.T @ qty
53
54 n_pens = len(S_list)
55 sp_vals = [1.0] * n_pens
56 free_indices = []
57 if sp_fixed is not None and len(sp_fixed) == n_pens:
58 for idx, val in enumerate(sp_fixed):
59 if val is not None:
60 sp_vals[idx] = float(val)
61 else:
62 free_indices.append(idx)
63 else:
64 free_indices = list(range(n_pens))
65
66 pen_bounds: List[Tuple[int, int]] = []
67 pen_sub: List[np.ndarray] = []
68 rank_S_list: List[int] = []
69 log_det_const: List[float] = []
70 supports: List[np.ndarray] = []
71 for S in S_list:
72 nz = np.where(np.diag(S) > 0)[0]
73 s0, s1 = (int(nz[0]), int(nz[-1]) + 1) if len(nz) > 0 else (0, 0)
74 pen_bounds.append((s0, s1))
75 sub = S[s0:s1, s0:s1]
76 pen_sub.append(sub)
77 supports.append(nz)
78 if s1 > s0:
79 ev = linalg.eigvalsh(sub, check_finite=False)
80 pos = ev[ev > max(float(ev[-1]) * 1e-10, 1e-14)]
81 rank_S_list.append(len(pos))
82 log_det_const.append(
83 float(np.sum(np.log(pos))) if len(pos) > 0 else 0.0
84 )
85 else:
86 rank_S_list.append(0)
87 log_det_const.append(0.0)
88
89 overlapping = [False] * n_pens
90 for i in range(n_pens):
91 for j in range(i + 1, n_pens):
92 if (
93 pen_bounds[i][0] < pen_bounds[j][1]
94 and pen_bounds[j][0] < pen_bounds[i][1]
95 and np.intersect1d(supports[i], supports[j]).size > 0
96 ):
97 overlapping[i] = overlapping[j] = True
98
99 overlap_groups: List[Tuple[List[int], int, Dict[int, np.ndarray]]] = []
100 visited = [False] * n_pens
101 for i in range(n_pens):
102 if overlapping[i] and not visited[i]:
103 grp = [i]
104 visited[i] = True
105 changed = True
106 while changed:
107 changed = False
108 for j in range(n_pens):
109 if overlapping[j] and not visited[j]:
110 if any(
111 np.intersect1d(supports[g], supports[j]).size > 0
112 for g in grp
113 ):
114 grp.append(j)
115 visited[j] = True
116 changed = True
117 g_s0 = min(pen_bounds[g][0] for g in grp)
118 g_s1 = max(pen_bounds[g][1] for g in grp)
119 S_sum = np.zeros((g_s1 - g_s0, g_s1 - g_s0), dtype=float)
120 for g_i in grp:
121 s0, s1 = pen_bounds[g_i]
122 S_sum[s0 - g_s0 : s1 - g_s0, s0 - g_s0 : s1 - g_s0] += pen_sub[
123 g_i
124 ]
125 ev0, U0 = linalg.eigh(S_sum, check_finite=False)
126 U_pos = U0[:, ev0 > max(float(ev0[-1]) * 1e-10, 1e-14)]
127 proj_subs = {
128 g_i: U_pos[
129 pen_bounds[g_i][0] - g_s0 : pen_bounds[g_i][1] - g_s0, :
130 ].T
131 @ pen_sub[g_i]
132 @ U_pos[
133 pen_bounds[g_i][0] - g_s0 : pen_bounds[g_i][1] - g_s0, :
134 ]
135 for g_i in grp
136 }
137 overlap_groups.append((grp, U_pos.shape[1], proj_subs))
138
139 total_pen_rank = sum(
140 rank_S_list[i] for i in range(n_pens) if not overlapping[i]
141 ) + sum(r_dim for _, r_dim, _ in overlap_groups)
142 M_p = max(p - total_pen_rank, 0)
143
144 eye_p = np.eye(p, dtype=float)
145
146 def _eval_criterion(
147 rho_vec: np.ndarray,
148 compute_details: bool = False,
149 compute_grad: bool = False,
150 ) -> Any:
151 current_sp = list(sp_vals)
152 for i, idx in enumerate(free_indices):
153 current_sp[idx] = float(np.exp(rho_vec[i]))
154
155 S_total = np.zeros((p, p), dtype=float)
156 for (s0, s1), sub, l_val in zip(pen_bounds, pen_sub, current_sp):
157 if s1 > s0:
158 S_total[s0:s1, s0:s1] += l_val * sub
159
160 H = RtR + S_total
161 H.flat[:: p + 1] += 1e-12
162
163 try:
164 c, lower = linalg.cho_factor(H, check_finite=False)
165 coef = linalg.cho_solve((c, lower), Xty, check_finite=False)
166 log_det_H = 2.0 * float(np.sum(np.log(np.diag(c))))
167 except Exception:
168 coef = np.linalg.solve(H, Xty)
169 log_det_H = float(np.linalg.slogdet(H)[1])
170 c, lower = None, False
171
172 r_diff = qty - R @ coef
173 rss = max(0.0, ss_orth + float(r_diff @ r_diff))
174 Dp = rss + float(coef @ (S_total @ coef))
175
176 H_inv = None
177 if compute_details or compute_grad:
178 if c is not None:
179 try:
180 H_inv = linalg.cho_solve(
181 (c, lower), eye_p, check_finite=False
182 )
183 except Exception:
184 H_inv = np.linalg.pinv(H)
185 else:
186 H_inv = np.linalg.pinv(H)
187
188 edf = 0.0
189 grad = np.zeros(len(free_indices), dtype=float) if compute_grad else None
190
191 if method == "REML":
192 df_reml = max(n - M_p, 1)
193 scale_est = max(Dp / df_reml, 1e-12)
194 log_det_S = sum(
195 rank_S_list[i] * np.log(current_sp[i]) + log_det_const[i]
196 for i in range(n_pens)
197 if not overlapping[i]
198 )
199 S_inv_proj: Dict[int, np.ndarray] = {}
200 for grp, r_dim, proj_subs in overlap_groups:
201 S_tilde = sum(current_sp[g_i] * proj_subs[g_i] for g_i in grp)
202 d = np.sqrt(np.maximum(np.diag(S_tilde), 1e-300))
203 S_bal = S_tilde / (d[:, None] * d[None, :])
204 try:
205 c_s, low_s = linalg.cho_factor(S_bal, check_finite=False)
206 log_det_S += 2.0 * float(
207 np.sum(np.log(np.diag(c_s))) + np.sum(np.log(d))
208 )
209 if compute_grad:
210 S_bal_inv = linalg.cho_solve(
211 (c_s, low_s), np.eye(r_dim), check_finite=False
212 )
213 S_inv = S_bal_inv / (d[:, None] * d[None, :])
214 for g_i in grp:
215 S_inv_proj[g_i] = S_inv
216 except Exception:
217 ev_s, u_s = linalg.eigh(S_tilde, check_finite=False)
218 pos_s = np.maximum(ev_s, 1e-300)
219 log_det_S += float(np.sum(np.log(pos_s)))
220 if compute_grad:
221 S_inv = (u_s / pos_s[None, :]) @ u_s.T
222 for g_i in grp:
223 S_inv_proj[g_i] = S_inv
224
225 score = (
226 Dp / (2.0 * scale_est)
227 + 0.5 * df_reml * np.log(2.0 * np.pi * scale_est)
228 + 0.5 * log_det_H
229 - 0.5 * log_det_S
230 )
231 if compute_grad and H_inv is not None and grad is not None:
232 for k_i, idx in enumerate(free_indices):
233 s0, s1 = pen_bounds[idx]
234 sub = pen_sub[idx]
235 lam_i = current_sp[idx]
236 b_sub = coef[s0:s1]
237 dDp = lam_i * float(b_sub @ sub @ b_sub)
238 dlogH = lam_i * float(np.sum(H_inv[s0:s1, s0:s1] * sub))
239 if not overlapping[idx]:
240 dlogS = float(rank_S_list[idx])
241 else:
242 dlogS = lam_i * float(
243 np.sum(
244 S_inv_proj[idx]
245 * next(
246 ps[idx]
247 for grp, _, ps in overlap_groups
248 if idx in grp
249 )
250 )
251 )
252 grad[k_i] = 0.5 * (dDp / scale_est + dlogH - dlogS)
253 else:
254 if H_inv is not None:
255 edf = float(np.sum(H_inv * RtR))
256 elif c is not None:
257 try:
258 Z_Rt = linalg.solve_triangular(
259 c,
260 R.T,
261 lower=lower,
262 trans="N" if lower else "T",
263 check_finite=False,
264 )
265 edf = float(np.sum(Z_Rt**2))
266 except Exception:
267 H_inv = linalg.cho_solve(
268 (c, lower), eye_p, check_finite=False
269 )
270 edf = float(np.sum(H_inv * RtR))
271 else:
272 H_inv = np.linalg.pinv(H)
273 edf = float(np.sum(H_inv * RtR))
274
275 denom = max(n - gamma * edf, 1e-6)
276 score = (n * rss) / (denom**2)
277 scale_est = rss / max(n - edf, 1)
278 if compute_grad and H_inv is not None and grad is not None:
279 res_grad_vec = H_inv @ (RtR @ coef - Xty)
280 H_RtR_H = (H_inv @ RtR) @ H_inv
281 for k_i, idx in enumerate(free_indices):
282 s0, s1 = pen_bounds[idx]
283 sub = pen_sub[idx]
284 lam_i = current_sp[idx]
285 drss = (
286 -2.0
287 * lam_i
288 * float(res_grad_vec[s0:s1] @ (sub @ coef[s0:s1]))
289 )
290 dedf = -lam_i * float(np.sum(H_RtR_H[s0:s1, s0:s1] * sub))
291 grad[k_i] = (n / (denom**2)) * drss + (
292 2.0 * n * rss * gamma / (denom**3)
293 ) * dedf
294
295 if compute_details and H_inv is not None and edf == 0.0:
296 edf = float(np.sum(H_inv * RtR))
297
298 if compute_grad:
299 return score, grad
300 return score, coef, H_inv, rss, scale_est, edf
301
302 if len(free_indices) == 0:
303 best_score, best_coef, best_H_inv, best_rss, _, best_edf = (
304 _eval_criterion(np.zeros(0), compute_details=True)
305 )
306 best_scale = best_rss / max(n - best_edf, 1)
307 best_sp = sp_vals
308 else:
309 n_free = len(free_indices)
310 candidates = [np.full(n_free, v) for v in (-8.0, -4.0, 0.0, 4.0, 8.0)]
311 rho_init = min(candidates, key=lambda r: _eval_criterion(r)[0])
312 opt = optimize.minimize(
313 lambda r: _eval_criterion(r, compute_grad=True),
314 rho_init,
315 jac=True,
316 method="L-BFGS-B",
317 bounds=[(-15.0, 15.0)] * n_free,
318 options={"ftol": 1e-14, "gtol": 1e-10},
319 )
320 best_score, best_coef, best_H_inv, best_rss, best_scale, best_edf = (
321 _eval_criterion(opt.x, compute_details=True)
322 )
323 best_sp = list(sp_vals)
324 for i, idx in enumerate(free_indices):
325 best_sp[idx] = float(np.exp(opt.x[i]))
326
327 scale_val = best_scale
328 Vp = best_H_inv * scale_val
329 H_inv_RtR = best_H_inv @ RtR
330 Ve = (H_inv_RtR @ best_H_inv) * scale_val
331 se = np.sqrt(np.maximum(np.diag(Vp), 1e-14))
332
333 edf_per_coef = np.diag(H_inv_RtR)
334 total_edf = float(np.sum(edf_per_coef))
335
336 edf_per_term: List[float] = []
337 p_values_per_term: List[float] = []
338 f_stat_per_term: List[float] = []
339
340 if coef_slices is not None:
341 df2 = max(n - total_edf, 1.0)
342 for start, end in coef_slices:
343 term_edf = float(np.sum(edf_per_coef[start:end]))
344 edf_per_term.append(term_edf)
345 beta_sub = best_coef[start:end]
346 V_sub = Vp[start:end, start:end]
347 try:
348 stat = float(
349 (beta_sub.T @ linalg.pinv(V_sub) @ beta_sub)
350 / max(term_edf, 1e-6)
351 )
352 pval = float(1.0 - stats.f.cdf(stat, max(term_edf, 1.0), df2))
353 except Exception:
354 stat, pval = 0.0, 1.0
355 f_stat_per_term.append(stat)
356 p_values_per_term.append(pval)
357
358 dev_explained = max(0.0, 1.0 - best_rss / max(null_deviance, 1e-12))
359 r_sq_adj = max(
360 0.0,
361 1.0
362 - (best_rss / max(n - total_edf, 1))
363 / (null_deviance / max(n - 1, 1)),
364 )
365
366 gcv_score = (n * best_rss) / (max(n - gamma * total_edf, 1.0) ** 2)
367 loglik = (
368 -0.5 * n * np.log(2.0 * np.pi * scale_val)
369 - 0.5 * best_rss / scale_val
370 + 0.5 * log_w_sum
371 )
372
373 return {
374 "coef": best_coef,
375 "sp": best_sp,
376 "scale": scale_val,
377 "Vp": Vp,
378 "Ve": Ve,
379 "se": se,
380 "edof": total_edf,
381 "edof_per_coef": edf_per_coef,
382 "edf_per_term": edf_per_term,
383 "f_stat_per_term": f_stat_per_term,
384 "p_values_per_term": p_values_per_term,
385 "rss": best_rss,
386 "deviance": best_rss,
387 "null_deviance": null_deviance,
388 "deviance_explained": dev_explained,
389 "r_sq_adj": r_sq_adj,
390 "GCV": gcv_score,
391 "score": best_score,
392 "method": method,
393 "loglik": loglik,
394 "log_weights_sum": log_w_sum,
395 "aic": -2.0 * loglik + 2.0 * (total_edf + 1),
396 "bic": -2.0 * loglik + np.log(n) * (total_edf + 1),
397 "n_samples": n,
398 "H_inv": best_H_inv,
399 "RtR": RtR,
400 }
magli.wood._terms¶
Term classes for Wood GAM backend: spline, tensor product, factor, and linear terms.
- magli.wood._terms.instantiate_wood_terms(configs: List[Dict[str, Any]]) List[_Term][source]¶
Instantiate Wood term objects from parsed formula config dicts.
1"""Term classes for Wood GAM backend: spline, tensor product, factor, and linear terms."""
2
3from __future__ import annotations
4
5from typing import Any, Dict, List, Optional, Tuple, Union
6
7import numpy as np
8
9from .._base import _khatri_rao
10from ..eilers._terms import _FactorTerm, _LinearTerm, _Term as _WoodTerm
11from ._smooths import (
12 absorb_sum_to_zero_constraint,
13 build_smooth_basis,
14 predict_smooth_basis,
15 rescale_penalty,
16)
17
18
19class _WoodSplineTerm(_WoodTerm):
20 """Univariate smooth term s(x, bs='tp', k=10)."""
21
22 def __init__(
23 self,
24 feature: int,
25 k: int = 10,
26 bs: str = "tp",
27 sp: Optional[float] = None,
28 m: Optional[Union[int, Tuple[int, int]]] = None,
29 knots: Optional[np.ndarray] = None,
30 constraint: Optional[str] = None,
31 ) -> None:
32 self.feature = int(feature)
33 self.k = int(k)
34 self.bs = bs.lower()
35 if constraint == "periodic" and self.bs == "tp":
36 self.bs = "cc"
37 self.sp = [float(sp)] if sp is not None else None
38 self.m = m
39 self.knots = knots
40 self._Z: Optional[np.ndarray] = None
41 self._penalties: List[np.ndarray] = []
42 self._meta: Dict[str, Any] = {}
43 self._cached_col: Optional[np.ndarray] = None
44 self._cached_X: Optional[np.ndarray] = None
45 self.n_coefs = 0
46
47 def compile(self, x: np.ndarray) -> None:
48 col = np.asarray(x[:, self.feature], dtype=float).ravel()
49 X_raw, S_list, self._meta = build_smooth_basis(
50 col, bs=self.bs, k=self.k, knots=self.knots, m=self.m
51 )
52 S_scaled, _ = rescale_penalty(X_raw, S_list)
53 X_abs, self._penalties, self._Z = absorb_sum_to_zero_constraint(
54 X_raw, S_scaled
55 )
56 self.n_coefs = X_abs.shape[1]
57 self._cached_col = col
58 self._cached_X = X_abs
59
60 def build_columns(self, x: np.ndarray) -> np.ndarray:
61 col = np.asarray(x[:, self.feature], dtype=float).ravel()
62 if self._Z is None:
63 self.compile(x)
64 if (
65 self._cached_col is not None
66 and col.shape == self._cached_col.shape
67 and np.array_equal(col, self._cached_col)
68 ):
69 return self._cached_X # type: ignore[return-value]
70 return predict_smooth_basis(col, self._meta) @ self._Z
71
72 def build_penalties(self) -> List[np.ndarray]:
73 return [p.copy() for p in self._penalties]
74
75
76class _WoodTensorTerm(_WoodTerm):
77 """Tensor product smooth te(x1, x2, bs='cr', k=5)."""
78
79 def __init__(
80 self,
81 features: List[int],
82 k: Union[int, List[int]] = 5,
83 bs: Union[str, List[str]] = "cr",
84 sp: Optional[Union[float, List[float]]] = None,
85 ) -> None:
86 self.features = [int(f) for f in features]
87 self.d = len(self.features)
88 if self.d < 2:
89 raise ValueError(
90 "Tensor product term requires at least 2 features."
91 )
92
93 self.k_per = (
94 [k] * self.d if isinstance(k, int) else [int(v) for v in k]
95 )
96 self.bs_per = (
97 [bs] * self.d if isinstance(bs, str) else [str(v) for v in bs]
98 )
99 self.sp = (
100 None
101 if sp is None
102 else [float(sp)] * self.d
103 if isinstance(sp, (int, float))
104 else [float(v) for v in sp]
105 )
106 self._Z: Optional[np.ndarray] = None
107 self._penalties: List[np.ndarray] = []
108 self._marginal_meta: List[Dict] = []
109 self._cached_cols: Optional[np.ndarray] = None
110 self._cached_X: Optional[np.ndarray] = None
111 self.n_coefs = 0
112
113 def compile(self, x: np.ndarray) -> None:
114 marg_X, marg_S, self._marginal_meta = [], [], []
115 cols = np.asarray(x[:, self.features], dtype=float)
116 for i in range(self.d):
117 col = cols[:, i].ravel()
118 Xi, Si, meta_i = build_smooth_basis(
119 col, bs=self.bs_per[i], k=self.k_per[i]
120 )
121 ev_max = np.linalg.eigvalsh(Si[0])[-1]
122 if ev_max > 0:
123 Si[0] = Si[0] / ev_max
124 marg_X.append(Xi)
125 marg_S.append(Si[0])
126 self._marginal_meta.append(meta_i)
127
128 X_te = _khatri_rao(marg_X)
129 sizes = [m.shape[1] for m in marg_X]
130 TS_list = []
131 for i in range(self.d):
132 I_pre = np.eye(int(np.prod(sizes[:i])))
133 I_post = np.eye(int(np.prod(sizes[i + 1 :])))
134 P_i = np.kron(np.kron(I_pre, marg_S[i]), I_post)
135 TS_list.append(0.5 * (P_i + P_i.T))
136
137 S_scaled, _ = rescale_penalty(X_te, TS_list)
138 X_abs, self._penalties, self._Z = absorb_sum_to_zero_constraint(
139 X_te, S_scaled
140 )
141 self.n_coefs = X_abs.shape[1]
142 self._cached_cols = cols
143 self._cached_X = X_abs
144
145 def build_columns(self, x: np.ndarray) -> np.ndarray:
146 if self._Z is None:
147 self.compile(x)
148 cols = np.asarray(x[:, self.features], dtype=float)
149 if (
150 self._cached_cols is not None
151 and cols.shape == self._cached_cols.shape
152 and np.array_equal(cols, self._cached_cols)
153 ):
154 return self._cached_X # type: ignore[return-value]
155 marg_X = [
156 predict_smooth_basis(
157 cols[:, i].ravel(),
158 self._marginal_meta[i],
159 )
160 for i in range(self.d)
161 ]
162 return _khatri_rao(marg_X) @ self._Z
163
164 def build_penalties(self) -> List[np.ndarray]:
165 return [p.copy() for p in self._penalties]
166
167
168class _WoodFactorTerm(_FactorTerm):
169 """Categorical factor term f(feature)."""
170
171 def __init__(self, feature: int, coding: str = "dummy") -> None:
172 super().__init__(feature, lam=0.0, coding=coding)
173
174
175_WoodLinearTerm = _LinearTerm
176
177
178def instantiate_wood_terms(configs: List[Dict[str, Any]]) -> List[_WoodTerm]:
179 """Instantiate Wood term objects from parsed formula config dicts."""
180 terms: List[_WoodTerm] = []
181 for cfg in configs:
182 kind = cfg["type"]
183 feats = cfg["features"]
184 kw = dict(cfg["kwargs"])
185
186 k_val = kw.pop("k", kw.pop("n_splines", 10))
187 sp_val = kw.pop("sp", kw.pop("lam", None))
188 bs_val = kw.pop("bs", "tp")
189
190 if kind == "s":
191 terms.append(
192 _WoodSplineTerm(
193 feature=feats[0],
194 k=k_val,
195 bs=bs_val,
196 sp=sp_val,
197 m=kw.pop("m", None),
198 knots=kw.pop("knots", None),
199 constraint=kw.pop("constraint", None),
200 )
201 )
202 elif kind == "te":
203 te_k = k_val if isinstance(k_val, (int, list)) else 5
204 terms.append(
205 _WoodTensorTerm(
206 features=feats, k=te_k, bs=bs_val, sp=sp_val
207 )
208 )
209 elif kind == "f":
210 terms.append(
211 _WoodFactorTerm(
212 feature=feats[0], coding=kw.pop("coding", "dummy")
213 )
214 )
215 elif kind == "l":
216 terms.append(_WoodLinearTerm(feature=feats[0]))
217 else:
218 raise ValueError(f"Unknown term type '{kind}' in formula.")
219 return terms
magli.wood._lingam¶
Single-term penalised GAM: LinGAM class (Wood backend).
- class magli.wood._lingam.LinGAM(n_splines: int = 10, lam: float | None = None, spline_order: int = 3, constraint: str | None = None, k: int | None = None, sp: float | None = None, bs: str = 'tp', method: str = 'REML', gamma: float = 1.0)[source]¶
Penalised GAM for a single smooth term y = beta_0 + f(x) + eps (Wood backend).
- fit(x: ArrayLike, y: ArrayLike, weights: ArrayLike | None = None) LinGAM[source]¶
Fit the GAM model using Wood REML or GCV criteria.
1"""Single-term penalised GAM: LinGAM class (Wood backend)."""
2
3from __future__ import annotations
4
5from typing import Optional, Tuple
6
7import numpy as np
8from numpy.typing import ArrayLike
9
10from ._gam import GAMCore
11from ._terms import _WoodSplineTerm
12
13
14class LinGAM(GAMCore):
15 """Penalised GAM for a single smooth term y = beta_0 + f(x) + eps (Wood backend)."""
16
17 def __init__(
18 self,
19 n_splines: int = 10,
20 lam: Optional[float] = None,
21 spline_order: int = 3,
22 constraint: Optional[str] = None,
23 k: Optional[int] = None,
24 sp: Optional[float] = None,
25 bs: str = "tp",
26 method: str = "REML",
27 gamma: float = 1.0,
28 ) -> None:
29 if k is not None:
30 n_splines = int(k)
31 if sp is not None:
32 lam = float(sp)
33
34 if constraint == "periodic":
35 bs = "cc"
36 elif constraint is not None:
37 raise ValueError(
38 f"Unsupported constraint '{constraint}' for Wood backend. "
39 "Wood backend natively supports 'periodic' (bs='cc'). For soft shape "
40 "constraints (monotonicity, convexity), use backend='eilers'."
41 )
42
43 if n_splines < 3:
44 raise ValueError("n_splines (k) must be at least 3.")
45 if lam is not None and lam < 0:
46 raise ValueError(
47 "Smoothing parameter 'lam' / 'sp' must be non-negative."
48 )
49
50 self.constraint = constraint
51 self.edge_knots_: Optional[np.ndarray] = None
52 self.knots_: Optional[np.ndarray] = None
53
54 bs_lower = bs.lower()
55 m_param = (
56 (int(spline_order), 2)
57 if bs_lower in ("bs", "b")
58 else (int(spline_order) - 1, 2)
59 if bs_lower in ("ps", "cp")
60 else None
61 )
62 lam_val = float(lam) if lam is not None else None
63 super().__init__(
64 f"s(0, k={int(n_splines)}, bs={bs_lower!r}, sp={lam_val!r}, m={m_param!r})",
65 method=method,
66 gamma=gamma,
67 fit_intercept=True,
68 spline_order=int(spline_order),
69 )
70
71 @property
72 def _term(self) -> _WoodSplineTerm:
73 return self._terms[0] # type: ignore[return-value]
74
75 @property
76 def n_splines(self) -> int:
77 return self._term.k
78
79 @n_splines.setter
80 def n_splines(self, val: int) -> None:
81 self._term.k = int(val)
82 self._term_configs[0]["kwargs"]["k"] = int(val)
83
84 k = n_splines
85
86 @property
87 def lam(self) -> Optional[float]:
88 return self._term.sp[0] if self._term.sp else None
89
90 @lam.setter
91 def lam(self, val: Optional[float]) -> None:
92 self._term.sp = [float(val)] if val is not None else None
93 self._term_configs[0]["kwargs"]["sp"] = (
94 float(val) if val is not None else None
95 )
96
97 sp = lam
98
99 @property
100 def bs(self) -> str:
101 return self._term.bs
102
103 @property
104 def n_coefs(self) -> int:
105 """Effective number of model coefficients (smooth + intercept)."""
106 return len(self.coef_) if self.coef_ is not None else self.n_splines
107
108 @property
109 def design_matrix_(self) -> Optional[np.ndarray]:
110 return self._design_matrix_
111
112 def __repr__(self) -> str:
113 status = "fitted" if self.is_fitted else "unfitted"
114 sp_str = f"sp={self.lam:.4g}" if self.lam is not None else "sp=auto"
115 cons_str = (
116 f", constraint='{self.constraint}'" if self.constraint else ""
117 )
118 return (
119 f"LinGAM(k={self.n_splines}, bs='{self.bs}', {sp_str}, "
120 f"method='{self.method}'{cons_str}) [{status}]"
121 )
122
123 def _diagnostics_title(self) -> str:
124 return f"LinGAM Diagnostics (k={self.n_splines}, bs='{self.bs}', sp={self.lam:.4f})"
125
126 def fit(
127 self,
128 x: ArrayLike,
129 y: ArrayLike,
130 weights: Optional[ArrayLike] = None,
131 ) -> "LinGAM":
132 """Fit the GAM model using Wood REML or GCV criteria."""
133 x_arr = np.asarray(x, dtype=float).ravel()
134 y_arr = np.asarray(y, dtype=float).ravel()
135 if not np.all(np.isfinite(x_arr)) or not np.all(np.isfinite(y_arr)):
136 raise ValueError("Input data contains NaN or infinity.")
137 self.edge_knots_ = np.array(
138 [float(np.min(x_arr)), float(np.max(x_arr))]
139 )
140 cfg_sp = self._term_configs[0]["kwargs"].get("sp")
141 self._term.sp = [float(cfg_sp)] if cfg_sp is not None else None
142 super().fit(x_arr.reshape(-1, 1), y_arr, weights=weights)
143 meta = self._term._meta
144 xp = meta.get("xp") if meta.get("xp") is not None else meta.get("knots")
145 if xp is None:
146 xp = meta.get("Xu")
147 self.knots_ = (
148 np.asarray(xp)
149 if xp is not None
150 else np.linspace(
151 self.edge_knots_[0], self.edge_knots_[1], self.n_splines
152 )
153 )
154 return self
155
156 def gridsearch(
157 self,
158 x: ArrayLike,
159 y: ArrayLike,
160 n_splines: Optional[ArrayLike] = None,
161 lam: Optional[ArrayLike] = None,
162 gamma: float = 1.0,
163 weights: Optional[ArrayLike] = None,
164 ) -> "LinGAM":
165 """Select optimal basis dimension k (and optionally lam) via REML or GCV."""
166 x_arr = np.asarray(x, dtype=float).ravel()
167 y_arr = np.asarray(y, dtype=float).ravel()
168 k_list = [
169 int(v)
170 for v in np.atleast_1d(
171 np.arange(5, 25, 2) if n_splines is None else n_splines
172 )
173 ]
174 lam_list = (
175 [None]
176 if lam is None
177 else [float(v) for v in np.atleast_1d(lam)]
178 )
179 crit_key = "GCV" if len(k_list) > 1 else "score"
180 best_score = np.inf
181 best_cand: Optional[LinGAM] = None
182
183 for k_cand in k_list:
184 for l_cand in lam_list:
185 try:
186 cand = LinGAM(
187 n_splines=k_cand,
188 lam=l_cand,
189 bs=self.bs,
190 method=self.method,
191 gamma=gamma,
192 spline_order=self.spline_order,
193 ).fit(x_arr, y_arr, weights=weights)
194 if cand.statistics_[crit_key] < best_score:
195 best_score = cand.statistics_[crit_key]
196 best_cand = cand
197 except Exception:
198 continue
199
200 if best_cand is not None:
201 self._term_configs = best_cand._term_configs
202 self._terms = best_cand._terms
203 self._coef_slices = best_cand._coef_slices
204 self.coef_ = best_cand.coef_
205 self.statistics_ = best_cand.statistics_
206 self._design_matrix_ = best_cand._design_matrix_
207 self.edge_knots_ = best_cand.edge_knots_
208 self.knots_ = best_cand.knots_
209 return self
210
211 return self.fit(x_arr, y_arr, weights=weights)
212
213 def plot_decomposition(
214 self,
215 x: ArrayLike,
216 y: ArrayLike,
217 n_grid: int = 300,
218 figsize: Tuple[int, int] = (16, 5),
219 show_weighted: bool = False,
220 ) -> Tuple:
221 """Plot basis decomposition, coefficients, and overall fit."""
222 k_smooth = self._term.n_coefs
223 fit_title = (
224 "Weighted Basis Functions + Fit"
225 if show_weighted
226 else "Overall Model Fit"
227 )
228 return self._plot_lingam_decomposition(
229 x,
230 y,
231 k_smooth=k_smooth,
232 prefix="b",
233 basis_title=f"Basis Functions (dim={k_smooth})",
234 fit_title=fit_title,
235 suptitle=f"LinGAM Decomposition (k={self.n_splines}, bs='{self.bs}', sp={self.lam:.4f}, method={self.method})",
236 knots=self.knots_,
237 show_weighted=show_weighted,
238 n_grid=n_grid,
239 figsize=figsize,
240 )
magli.wood._gam¶
Multi-term GAM with formula interface and REML/GCV optimization (Wood backend).
- class magli.wood._gam.GAMCore(formula: str, method: str = 'REML', gamma: float = 1.0, fit_intercept: bool = True, spline_order: int = 3)[source]¶
Multi-term Generalized Additive Model with formula interface (Wood backend).
- fit(x: ArrayLike, y: ArrayLike, weights: ArrayLike | None = None) GAMCore[source]¶
Fit multi-term GAM model using REML or GCV criteria.
1"""Multi-term GAM with formula interface and REML/GCV optimization (Wood backend)."""
2
3from __future__ import annotations
4
5from itertools import product as cartesian_product
6from typing import Any, Dict, List, Optional, Tuple
7
8import numpy as np
9from numpy.typing import ArrayLike
10from scipy import stats
11
12from .._base import _GAMBase, _validate_weights
13from .._formula import parse_formula
14from ._optimizer import fit_wood_engine
15from ._terms import (
16 _WoodLinearTerm,
17 _WoodSplineTerm,
18 _WoodTensorTerm,
19 _WoodTerm,
20 instantiate_wood_terms,
21)
22
23
24class GAMCore(_GAMBase):
25 """Multi-term Generalized Additive Model with formula interface (Wood backend)."""
26
27 _scale_is_var = True
28
29 def __init__(
30 self,
31 formula: str,
32 method: str = "REML",
33 gamma: float = 1.0,
34 fit_intercept: bool = True,
35 spline_order: int = 3,
36 ) -> None:
37 self.formula = formula
38 self.method = method.upper()
39 self.gamma = float(gamma)
40 self.fit_intercept = bool(fit_intercept)
41 self.spline_order = int(spline_order)
42
43 configs = parse_formula(formula)
44 self._term_configs = configs
45 self._terms: List[_WoodTerm] = instantiate_wood_terms(configs)
46
47 self._coef_slices: List[Tuple[int, int]] = []
48 self.coef_: Optional[np.ndarray] = None
49 self.statistics_: Optional[Dict[str, Any]] = None
50 self._design_matrix_: Optional[np.ndarray] = None
51
52 def __repr__(self) -> str:
53 status = "fitted" if self.is_fitted else "unfitted"
54 return f"GAMCore(formula='{self.formula}', method='{self.method}') [{status}]"
55
56 def fit(
57 self,
58 x: ArrayLike,
59 y: ArrayLike,
60 weights: Optional[ArrayLike] = None,
61 ) -> "GAMCore":
62 """Fit multi-term GAM model using REML or GCV criteria."""
63 x_arr = np.asarray(x, dtype=float)
64 y_arr = np.asarray(y, dtype=float).ravel()
65 if not np.all(np.isfinite(x_arr)) or not np.all(np.isfinite(y_arr)):
66 raise ValueError("Input data contains NaN or infinity.")
67 if x_arr.ndim == 1:
68 x_arr = x_arr.reshape(-1, 1)
69 w_arr = _validate_weights(weights, len(y_arr))
70
71 self._compile(x_arr)
72 self._compute_coef_slices()
73 X = self._build_model_matrix(x_arr)
74 self._design_matrix_ = X
75 p = X.shape[1]
76
77 S_list: List[np.ndarray] = []
78 sp_fixed: List[Optional[float]] = []
79 pen_counts: List[int] = []
80
81 for i, term in enumerate(self._terms):
82 start, end = self._coef_slices[i]
83 pens = term.build_penalties()
84 pen_counts.append(len(pens))
85 term_sps = term.sp if term.sp is not None else [None] * len(pens)
86 for pen_idx, P_sub in enumerate(pens):
87 P_full = np.zeros((p, p), dtype=float)
88 P_full[start:end, start:end] = P_sub
89 S_list.append(P_full)
90 sp_fixed.append(
91 term_sps[pen_idx] if pen_idx < len(term_sps) else None
92 )
93
94 fit_res = fit_wood_engine(
95 X=X,
96 y=y_arr,
97 S_list=S_list,
98 sp_fixed=(
99 sp_fixed if any(sp is not None for sp in sp_fixed) else None
100 ),
101 method=self.method,
102 gamma=self.gamma,
103 coef_slices=self._coef_slices,
104 M_p=1 if self.fit_intercept else 0,
105 weights=w_arr,
106 )
107
108 self.coef_ = fit_res["coef"]
109
110 sp_idx = 0
111 for term, n_pens in zip(self._terms, pen_counts):
112 if n_pens > 0:
113 term.sp = fit_res["sp"][sp_idx : sp_idx + n_pens]
114 sp_idx += n_pens
115
116 self.statistics_ = {
117 "edof": fit_res["edof"],
118 "edof_per_coef": fit_res["edof_per_coef"],
119 "scale": fit_res["scale"],
120 "cov": fit_res["Vp"],
121 "Vp": fit_res["Vp"],
122 "Ve": fit_res["Ve"],
123 "se": fit_res["se"],
124 "rss": fit_res["rss"],
125 "deviance": fit_res["deviance"],
126 "null_deviance": fit_res["null_deviance"],
127 "deviance_explained": fit_res["deviance_explained"],
128 "r_sq_adj": fit_res["r_sq_adj"],
129 "n_samples": fit_res["n_samples"],
130 "GCV": fit_res["GCV"],
131 "score": fit_res["score"],
132 "method": fit_res["method"],
133 "sp": fit_res["sp"],
134 "edf_per_term": fit_res["edf_per_term"],
135 "f_stat_per_term": fit_res["f_stat_per_term"],
136 "p_values_per_term": fit_res["p_values_per_term"],
137 "log_weights_sum": fit_res["log_weights_sum"],
138 }
139 return self
140
141 def gridsearch(
142 self,
143 x: ArrayLike,
144 y: ArrayLike,
145 n_splines_grids: Optional[Any] = None,
146 gamma: float = 1.0,
147 n_splines: Optional[Any] = None,
148 search_terms: Optional[Any] = None,
149 weights: Optional[ArrayLike] = None,
150 **kwargs: Any,
151 ) -> "GAMCore":
152 """Search over candidate basis dimensions k for smooth terms via GCV."""
153 x_arr = np.asarray(x, dtype=float)
154 y_arr = np.asarray(y, dtype=float).ravel()
155 if x_arr.ndim == 1:
156 x_arr = x_arr.reshape(-1, 1)
157 w_arr = _validate_weights(weights, len(y_arr))
158
159 if n_splines_grids is None:
160 n_splines_grids = n_splines if n_splines is not None else kwargs.get("k")
161
162 all_smooth = [
163 i
164 for i, t in enumerate(self._terms)
165 if isinstance(t, (_WoodSplineTerm, _WoodTensorTerm))
166 ]
167 if not all_smooth:
168 return self.fit(x_arr, y_arr, weights=w_arr)
169
170 if search_terms is None:
171 smooth_indices = all_smooth
172 else:
173 items = [search_terms] if isinstance(search_terms, (int, str)) else list(search_terms)
174 allowed = set()
175 for item in items:
176 if isinstance(item, int):
177 allowed.add(item)
178 elif item == "all":
179 allowed.update(all_smooth)
180 elif item == "s":
181 allowed.update(i for i, t in enumerate(self._terms) if isinstance(t, _WoodSplineTerm))
182 elif item == "te":
183 allowed.update(i for i, t in enumerate(self._terms) if isinstance(t, _WoodTensorTerm))
184 smooth_indices = [i for i in all_smooth if i in allowed]
185 if not smooth_indices:
186 return self.fit(x_arr, y_arr, weights=w_arr)
187
188 if n_splines_grids is None:
189 grids = [
190 np.array([4, 6, 8])
191 if isinstance(self._terms[i], _WoodTensorTerm)
192 else np.arange(6, 16, 2)
193 for i in smooth_indices
194 ]
195 elif isinstance(n_splines_grids, list) and len(n_splines_grids) > 0 and isinstance(n_splines_grids[0], (list, tuple, np.ndarray)):
196 grids = [np.atleast_1d(g) for g in n_splines_grids]
197 else:
198 grids = [np.atleast_1d(n_splines_grids)] * len(smooth_indices)
199
200 best_gcv = np.inf
201 best_model: Optional[GAMCore] = None
202
203 for combo in cartesian_product(*grids):
204 configs_copy = [dict(cfg) for cfg in self._term_configs]
205 for idx_in_smooth, term_i in enumerate(smooth_indices):
206 cfg = configs_copy[term_i]
207 cfg["kwargs"] = {
208 **cfg.get("kwargs", {}),
209 "k": int(combo[idx_in_smooth]),
210 }
211 cand_model = GAMCore(
212 formula=self.formula,
213 method=self.method,
214 gamma=gamma,
215 fit_intercept=self.fit_intercept,
216 spline_order=self.spline_order,
217 )
218 cand_model._terms = instantiate_wood_terms(configs_copy)
219 cand_model._term_configs = configs_copy
220 try:
221 cand_model.fit(x_arr, y_arr, weights=w_arr)
222 if cand_model.statistics_["GCV"] < best_gcv:
223 best_gcv = cand_model.statistics_["GCV"]
224 best_model = cand_model
225 except Exception:
226 continue
227
228 if best_model is not None:
229 self._term_configs = best_model._term_configs
230 self._terms = best_model._terms
231 self._coef_slices = best_model._coef_slices
232 self.coef_ = best_model.coef_
233 self.statistics_ = best_model.statistics_
234 self._design_matrix_ = best_model._design_matrix_
235 return self
236
237 return self.fit(x_arr, y_arr, weights=w_arr)
238
239 def summary(self) -> str:
240 """Formatted summary table for the fitted Wood GAM."""
241 self._check_fitted()
242 n = self.statistics_["n_samples"]
243 score_name = "REML" if self.method == "REML" else "GCV"
244 score_val = self.statistics_["score"]
245
246 lines = [
247 "",
248 "Family: gaussian ",
249 "Link function: identity ",
250 "",
251 f"Formula: {self.formula}",
252 "",
253 ]
254
255 param_rows = []
256 df_res = max(n - self.statistics_["edof"], 1.0)
257 if self.fit_intercept:
258 b0 = self.coef_[-1]
259 se0 = self.statistics_["se"][-1]
260 t0 = b0 / max(se0, 1e-12)
261 p0 = 2.0 * (1.0 - stats.t.cdf(np.abs(t0), df=df_res))
262 param_rows.append(("(Intercept)", b0, se0, t0, p0))
263
264 for i, t in enumerate(self._terms):
265 if isinstance(t, _WoodLinearTerm):
266 start, _ = self._coef_slices[i]
267 b = self.coef_[start]
268 se_b = self.statistics_["se"][start]
269 tb = b / max(se_b, 1e-12)
270 pb = 2.0 * (1.0 - stats.t.cdf(np.abs(tb), df=df_res))
271 param_rows.append((f"l({t.feature})", b, se_b, tb, pb))
272
273 if param_rows:
274 lines.append("Parametric coefficients:")
275 lines.append(" Estimate Std. Error t value Pr(>|t|)")
276 for name, b, se_b, tb, pb in param_rows:
277 lines.append(
278 f"{name:11s} {b:10.4f} {se_b:10.4f} {tb:7.3f} {pb:8.4g}"
279 )
280 lines.append("")
281
282 smooth_rows = []
283 for i, t in enumerate(self._terms):
284 if isinstance(t, (_WoodSplineTerm, _WoodTensorTerm)):
285 if isinstance(t, _WoodSplineTerm):
286 label, ref_df = f"s({t.feature})", t.k - 1
287 else:
288 label = f"te({','.join(str(f) for f in t.features)})"
289 ref_df = t.n_coefs
290 smooth_rows.append(
291 (
292 label,
293 self.statistics_["edf_per_term"][i],
294 ref_df,
295 self.statistics_["f_stat_per_term"][i],
296 self.statistics_["p_values_per_term"][i],
297 )
298 )
299
300 if smooth_rows:
301 lines.append("Approximate significance of smooth terms:")
302 lines.append(" edf Ref.df F p-value")
303 for label, edf, ref_df, f_stat, pval in smooth_rows:
304 lines.append(
305 f"{label:7s} {edf:6.2f} {ref_df:6d} {f_stat:6.3f} {pval:8.4g}"
306 )
307 lines.append("")
308
309 lines.extend(
310 [
311 f"R-sq.(adj) = {self.statistics_['r_sq_adj']:.4f} Deviance explained = {self.statistics_['deviance_explained'] * 100:.1f}%",
312 f"{score_name} = {score_val:10.4f} Scale est. = {self.statistics_['scale']:.4f} n = {n}",
313 "",
314 ]
315 )
316 return "\n".join(lines)
317
318 def plot_decomposition(
319 self,
320 x: ArrayLike,
321 y: ArrayLike,
322 n_grid: int = 200,
323 figsize: Optional[Tuple[int, int]] = None,
324 feature_names: Optional[List[str]] = None,
325 ) -> Tuple:
326 """Plot partial dependence for all terms in the GAM."""
327 return self._plot_gamcore_decomposition(
328 x,
329 y,
330 n_grid=n_grid,
331 extrap_frac=0.0,
332 figsize=figsize,
333 feature_names=feature_names,
334 )