"""
Spline and other smoother classes for Generalized Additive Models
Author: Luca Puggini
Author: Josef Perktold
Created on Fri Jun 5 16:32:00 2015
"""
from statsmodels.compat.numpy import inplace_reshape
from statsmodels.compat.patsy import get_all_sorted_knots
from statsmodels.compat.python import with_metaclass
# import useful only for development
from abc import ABCMeta, abstractmethod
import numpy as np
import pandas as pd
from statsmodels.formula._manager import FormulaManager
from statsmodels.tools.linalg import transf_constraints
# Obtain b splines from patsy
def _equally_spaced_knots(x, df):
"""
Compute equally spaced knots over the range of x
Parameters
----------
x : ndarray
1-D array of the underlying explanatory variable.
df : int
Number of knots to compute is ``df - 2``.
Returns
-------
knots : ndarray
1-D array of equally spaced knots between the min and max of `x`.
"""
n_knots = df - 2
x_min = x.min()
x_max = x.max()
knots = np.linspace(x_min, x_max, n_knots)
return knots
def _R_compat_quantile(x, probs):
"""
Compute quantiles of x in a way that is compatible with R's quantile
Parameters
----------
x : ndarray
Data for which to compute quantiles.
probs : array_like
Probabilities at which to compute quantiles, in [0, 1].
Returns
-------
quantiles : ndarray
Quantiles of `x` at `probs`, with the same shape as `probs`.
"""
# return np.percentile(x, 100 * np.asarray(probs))
probs = np.asarray(probs)
quantiles = np.asarray([np.percentile(x, 100 * prob)
for prob in probs.ravel(order="C")])
return quantiles.reshape(probs.shape, order="C")
# FIXME: is this copy/pasted? If so, why do we need it? If not, get
# rid of the try/except for scipy import
# from patsy splines.py
def _eval_bspline_basis(x, knots, degree, deriv="all", include_intercept=True):
"""
Evaluate a B-spline basis, and optionally its derivatives, at x
Parameters
----------
x : array_like
Points at which to evaluate the basis.
knots : array_like
1-D array of knots. Assumed to already be pre-processed, e.g.
duplicate copies of the boundary knots should already be included.
degree : int
Degree of the spline.
deriv : {"all", 0, 1, 2}
Which derivative of the basis to return. If "all", a tuple with
the basis and its first and second derivatives is returned.
include_intercept : bool
If False, the first basis column (the constant) is dropped.
Returns
-------
basis : ndarray or tuple of ndarray
If `deriv` is "all", a tuple ``(basis, der1_basis, der2_basis)``
of design matrices for the spline basis, its first, and its
second derivative. Otherwise, the single design matrix for the
requested derivative.
"""
try:
from scipy.interpolate import splev
except ImportError as err:
raise ImportError("spline functionality requires scipy") from err
# 'knots' are assumed to be already pre-processed. E.g. usually you
# want to include duplicate copies of boundary knots; you should do
# that *before* calling this constructor.
knots = np.sort(np.atleast_1d(np.asarray(knots, dtype=float)))
assert knots.ndim == 1
degree = int(degree)
x = np.atleast_1d(x)
if x.ndim == 2 and x.shape[1] == 1:
x = x[:, 0]
assert x.ndim == 1
# XX FIXME: when points fall outside of the boundaries, splev and R seem
# to handle them differently. I do not know why yet. So until we understand
# this and decide what to do with it, I'm going to play it safe and
# disallow such points.
if np.min(x) < np.min(knots) or np.max(x) > np.max(knots):
raise NotImplementedError("some data points fall outside the "
"outermost knots, and I'm not sure how "
"to handle them. (Patches accepted!)")
# Thanks to Charles Harris for explaining splev. It's not well
# documented, but basically it computes an arbitrary b-spline basis
# given knots and degree on some specificed points (or derivatives
# thereof, but we do not use that functionality), and then returns some
# linear combination of these basis functions. To get out the basis
# functions themselves, we use linear combinations like [1, 0, 0], [0,
# 1, 0], [0, 0, 1].
# NB: This probably makes it rather inefficient (though I have not checked
# to be sure -- maybe the fortran code actually skips computing the basis
# function for coefficients that are zero).
# Note: the order of a spline is the same as its degree + 1.
# Note: there are (len(knots) - order) basis functions.
k_const = 1 - int(include_intercept)
n_bases = len(knots) - (degree + 1) - k_const
if deriv in ["all", 0]:
basis = np.empty((x.shape[0], n_bases), dtype=float)
ret = basis
if deriv in ["all", 1]:
der1_basis = np.empty((x.shape[0], n_bases), dtype=float)
ret = der1_basis
if deriv in ["all", 2]:
der2_basis = np.empty((x.shape[0], n_bases), dtype=float)
ret = der2_basis
for i in range(n_bases):
coefs = np.zeros((n_bases + k_const,))
# we are skipping the first column of the basis to drop constant
coefs[i + k_const] = 1
ii = i
if deriv in ["all", 0]:
basis[:, ii] = splev(x, (knots, coefs, degree))
if deriv in ["all", 1]:
der1_basis[:, ii] = splev(x, (knots, coefs, degree), der=1)
if deriv in ["all", 2]:
der2_basis[:, ii] = splev(x, (knots, coefs, degree), der=2)
if deriv == "all":
return basis, der1_basis, der2_basis
else:
return ret
def compute_all_knots(x, df, degree):
"""
Compute the interior and boundary knots for a B-spline basis
Parameters
----------
x : array_like
Underlying explanatory variable used to place the interior knots.
df : int
Number of basis functions or degrees of freedom.
degree : int
Degree of the spline.
Returns
-------
all_knots : ndarray
1-D array of all knots, boundary knots repeated ``degree + 1``
times followed by the interior knots.
lower_bound : float
Lower boundary of `x`.
upper_bound : float
Upper boundary of `x`.
inner_knots : ndarray
1-D array of interior knots.
"""
order = degree + 1
n_inner_knots = df - order
lower_bound = np.min(x)
upper_bound = np.max(x)
knot_quantiles = np.linspace(0, 1, n_inner_knots + 2)[1:-1]
inner_knots = _R_compat_quantile(x, knot_quantiles)
all_knots = np.concatenate(([lower_bound, upper_bound] * order,
inner_knots))
return all_knots, lower_bound, upper_bound, inner_knots
def make_bsplines_basis(x, df, degree):
"""Make a spline basis for x."""
all_knots, _, _, _ = compute_all_knots(x, df, degree)
basis, der_basis, der2_basis = _eval_bspline_basis(x, all_knots, degree)
return basis, der_basis, der2_basis
def get_knots_bsplines(x=None, df=None, knots=None, degree=3,
spacing="quantile", lower_bound=None,
upper_bound=None, all_knots=None):
"""
Compute knots for use in B-splines
There are two main options for the knot placement
- quantile spacing with multiplicity of boundary knots
- equal spacing extended to boundary or exterior knots
The first corresponds to splines as used by patsy, the second is the
knot spacing for P-Splines.
Parameters
----------
x : ndarray, optional
Underlying explanatory variable used to compute quantile-based
interior knots. Required unless `all_knots` is provided.
df : int, optional
Number of basis functions or degrees of freedom, used together
with `degree` to determine the number of interior knots. Either
`df` or `knots` must be provided.
knots : array_like, optional
Interior knots. If provided, these are used directly instead of
being computed from `x`.
degree : int
Degree of the spline.
spacing : {"quantile", "equal"}
Method used to compute interior knots when `knots` is not
provided: quantiles of `x`, or equally spaced over its range.
lower_bound : float, optional
Location of the lower boundary knots. If None, uses the minimum
of `x`.
upper_bound : float, optional
Location of the upper boundary knots. If None, uses the maximum
of `x`.
all_knots : array_like, optional
If provided, these are returned unchanged and all other
arguments are ignored.
Returns
-------
all_knots : ndarray
1-D array of sorted knots, including boundary and interior knots.
"""
# based on patsy memorize_finish
if all_knots is not None:
return all_knots
x_min = x.min()
x_max = x.max()
if degree < 0:
raise ValueError(f"degree must be greater than 0 (not {degree!r})")
if int(degree) != degree:
raise ValueError(f"degree must be an integer (not {degree!r})")
# These are guaranteed to all be 1d vectors by the code above
# x = np.concatenate(tmp["xs"])
if df is None and knots is None:
raise ValueError("must specify either df or knots")
order = degree + 1
if df is not None:
n_inner_knots = df - order
if n_inner_knots < 0:
# We know that n_inner_knots is negative; if df were that much
# larger, it would have been zero, and things would work.
min_df = df - n_inner_knots
raise ValueError(
f"df={df!r} is too small for degree={degree!r}; must be >= {min_df}"
)
if knots is not None:
if len(knots) != n_inner_knots:
raise ValueError(f"df={df} with degree={degree!r} implies {n_inner_knots} knots, "
f"but {len(knots)} knots were provided")
elif spacing == "quantile":
# Need to compute inner knots
knot_quantiles = np.linspace(0, 1, n_inner_knots + 2)[1:-1]
inner_knots = _R_compat_quantile(x, knot_quantiles)
elif spacing == "equal":
# Need to compute inner knots
grid = np.linspace(0, 1, n_inner_knots + 2)[1:-1]
inner_knots = x_min + grid * (x_max - x_min)
diff_knots = inner_knots[1] - inner_knots[0]
else:
raise ValueError("incorrect option for spacing")
if knots is not None:
inner_knots = knots
if lower_bound is None:
lower_bound = np.min(x)
if upper_bound is None:
upper_bound = np.max(x)
if lower_bound > upper_bound:
raise ValueError(f"lower_bound > upper_bound ({lower_bound!r} > {upper_bound!r})")
inner_knots = np.asarray(inner_knots)
if inner_knots.ndim > 1:
raise ValueError("knots must be 1 dimensional")
if np.any(inner_knots < lower_bound):
raise ValueError(f"some knot values ({inner_knots[inner_knots < lower_bound]}) fall below lower bound "
f"({lower_bound!r})")
if np.any(inner_knots > upper_bound):
raise ValueError(f"some knot values ({inner_knots[inner_knots > upper_bound]}) fall above upper bound "
f"({upper_bound!r})")
if spacing == "equal":
diffs = np.arange(1, order + 1) * diff_knots
lower_knots = inner_knots[0] - diffs[::-1]
upper_knots = inner_knots[-1] + diffs
all_knots = np.concatenate((lower_knots, inner_knots, upper_knots))
else:
all_knots = np.concatenate(([lower_bound, upper_bound] * order,
inner_knots))
all_knots.sort()
return all_knots
def _get_integration_points(knots, k_points=3):
"""
Add points to each subinterval defined by knots
Inserts `k_points` between each two consecutive knots.
Parameters
----------
knots : array_like
1-D array of knots.
k_points : int
Number of points to insert between each pair of consecutive
knots.
Returns
-------
x : ndarray
1-D array of the unique knots together with the inserted points.
"""
k_points = k_points + 1
knots = np.unique(knots)
dxi = np.arange(k_points) / k_points
dxk = np.diff(knots)
dx = dxk[:, None] * dxi
x = np.concatenate(((knots[:-1, None] + dx).ravel(), [knots[-1]]))
return x
def get_covder2(smoother, k_points=3, integration_points=None,
skip_ctransf=False, deriv=2):
"""
Approximate integral of cross product of second derivative of smoother
This uses scipy.integrate simps to compute an approximation to the
integral of the smoother derivative cross-product at knots plus k_points
in between knots.
Parameters
----------
smoother : instance of univariate smoother class
Smoother for which the penalty matrix based on the second
derivative is computed. Needs a ``knots`` attribute and a
``transform`` method.
k_points : int
Number of points to insert between each pair of consecutive knots
when `integration_points` is not provided.
integration_points : array_like, optional
Points at which to evaluate the second derivative. If not
provided, points are generated from `smoother.knots` using
`k_points`.
skip_ctransf : bool
Whether to skip the constraint transform when evaluating the
smoother's derivative.
deriv : int
Derivative of the smoother basis used to compute the penalty
matrix.
Returns
-------
covd2 : ndarray
Approximation to the integral of the cross product of the
`deriv`-th derivative of the smoother basis functions.
"""
from scipy.integrate import simpson
knots = smoother.knots
if integration_points is None:
x = _get_integration_points(knots, k_points=k_points)
else:
x = integration_points
d2 = smoother.transform(x, deriv=deriv, skip_ctransf=skip_ctransf)
covd2 = simpson(d2[:, :, None] * d2[:, None, :], x=x, axis=0)
return covd2
# TODO: this function should be deleted
def make_poly_basis(x, degree, intercept=True):
"""
Given a vector x returns poly=(1, x, x^2, ..., x^degree)
and its first and second derivative.
Parameters
----------
x : ndarray
1-D array of the underlying explanatory variable.
degree : int
Highest power included in the polynomial basis.
intercept : bool
If False, the constant column is dropped from the basis and its
derivatives.
Returns
-------
basis : ndarray
Design matrix with columns ``x**i`` for ``i`` from 0 (or 1 if
`intercept` is False) to `degree`.
der_basis : ndarray
First derivative of `basis` with respect to `x`.
der2_basis : ndarray
Second derivative of `basis` with respect to `x`.
"""
if intercept:
start = 0
else:
start = 1
nobs = len(x)
basis = np.zeros(shape=(nobs, degree + 1 - start))
der_basis = np.zeros(shape=(nobs, degree + 1 - start))
der2_basis = np.zeros(shape=(nobs, degree + 1 - start))
for i in range(start, degree + 1):
basis[:, i - start] = x ** i
der_basis[:, i - start] = i * x ** (i - 1)
der2_basis[:, i - start] = i * (i - 1) * x ** (i - 2)
return basis, der_basis, der2_basis
# TODO: try to include other kinds of splines from patsy
# x = np.linspace(0, 1, 30)
# df = 10
# degree = 3
# from patsy.mgcv_cubic_splines import cc, cr, te
# all_knots, lower, upper, inner = compute_all_knots(x, df, degree)
# result = cc(x, df=df, knots=all_knots, lower_bound=lower, upper_bound=upper,
# constraints=None)
#
# import matplotlib.pyplot as plt
#
# result = np.array(result)
# print(result.shape)
# plt.plot(result.T)
# plt.show()
class UnivariateGamSmoother(with_metaclass(ABCMeta)):
"""
Base class for single smooth component
Parameters
----------
x : ndarray, 1-D
Underlying explanatory variable for the smooth term.
constraints : {None, str, array_like}
Constraints are used to transform the basis functions to satisfy
those constraints. ``constraints = 'center'`` applies a linear
transform to remove the constant and center the basis functions.
variable_name : str
The name for the underlying explanatory variable, x, used for
creating the column and parameter names for the basis functions.
"""
def __init__(self, x, constraints=None, variable_name="x"):
self.x = x
self.constraints = constraints
self.variable_name = variable_name
self.nobs, self.k_variables = len(x), 1
base4 = self._smooth_basis_for_single_variable()
if constraints == "center":
constraints = base4[0].mean(0)[None, :]
if constraints is not None and not isinstance(constraints, str):
ctransf = transf_constraints(constraints)
self.ctransf = ctransf
# subclasses might set ctransf directly
# only used if constraints is None
elif not hasattr(self, "ctransf"):
self.ctransf = None
self.basis, self.der_basis, self.der2_basis, self.cov_der2 = base4
if self.ctransf is not None:
ctransf = self.ctransf
# transform attributes that are not None
if base4[0] is not None:
self.basis = base4[0].dot(ctransf)
if base4[1] is not None:
self.der_basis = base4[1].dot(ctransf)
if base4[2] is not None:
self.der2_basis = base4[2].dot(ctransf)
if base4[3] is not None:
self.cov_der2 = ctransf.T.dot(base4[3]).dot(ctransf)
self.dim_basis = self.basis.shape[1]
self.col_names = [self.variable_name + "_s" + str(i)
for i in range(self.dim_basis)]
@abstractmethod
def _smooth_basis_for_single_variable(self):
return
class UnivariateGenericSmoother(UnivariateGamSmoother):
"""
Generic single smooth component
Parameters
----------
x : ndarray, 1-D
Underlying explanatory variable for the smooth term.
basis : ndarray
Design matrix of the basis functions.
der_basis : ndarray
First derivative of `basis`.
der2_basis : ndarray
Second derivative of `basis`.
cov_der2 : ndarray
Penalty matrix based on the second derivative.
variable_name : str
The name for the underlying explanatory variable, x, used for
creating the column and parameter names for the basis functions.
"""
def __init__(self, x, basis, der_basis, der2_basis, cov_der2,
variable_name="x"):
self.basis = basis
self.der_basis = der_basis
self.der2_basis = der2_basis
self.cov_der2 = cov_der2
super().__init__(
x, variable_name=variable_name)
def _smooth_basis_for_single_variable(self):
return self.basis, self.der_basis, self.der2_basis, self.cov_der2
class UnivariatePolynomialSmoother(UnivariateGamSmoother):
"""
Polynomial single smooth component
Parameters
----------
x : ndarray, 1-D
Underlying explanatory variable for the smooth term.
degree : int
Highest power included in the polynomial basis.
variable_name : str
The name for the underlying explanatory variable, x, used for
creating the column and parameter names for the basis functions.
"""
def __init__(self, x, degree, variable_name="x"):
self.degree = degree
super().__init__(
x, variable_name=variable_name)
def _smooth_basis_for_single_variable(self):
# TODO: unclear description
"""
Given a vector x returns poly=(1, x, x^2, ..., x^degree)
and its first and second derivative.
"""
basis = np.zeros(shape=(self.nobs, self.degree))
der_basis = np.zeros(shape=(self.nobs, self.degree))
der2_basis = np.zeros(shape=(self.nobs, self.degree))
for i in range(self.degree):
dg = i + 1
basis[:, i] = self.x ** dg
der_basis[:, i] = dg * self.x ** (dg - 1)
if dg > 1:
der2_basis[:, i] = dg * (dg - 1) * self.x ** (dg - 2)
else:
der2_basis[:, i] = 0
cov_der2 = np.dot(der2_basis.T, der2_basis)
return basis, der_basis, der2_basis, cov_der2
class UnivariateBSplines(UnivariateGamSmoother):
"""
B-Spline single smooth component
This creates and holds the B-Spline basis function for one
component.
Parameters
----------
x : ndarray, 1-D
underlying explanatory variable for smooth terms.
df : int
number of basis functions or degrees of freedom
degree : int
degree of the spline
include_intercept : bool
If False, then the basis functions are transformed so that they
do not include a constant. This avoids perfect collinearity if
a constant or several components are included in the model.
constraints : {None, str, array}
Constraints are used to transform the basis functions to satisfy
those constraints.
`constraints = 'center'` applies a linear transform to remove the
constant and center the basis functions.
variable_name : {None, str}
The name for the underlying explanatory variable, x, used in for
creating the column and parameter names for the basis functions.
covder2_kwds : {None, dict}
options for computing the penalty matrix from the second derivative
of the spline.
knot_kwds : {None, list[dict]}
option for the knot selection.
By default knots are selected in the same way as in patsy, however the
number of knots is independent of keeping or removing the constant.
Interior knot selection is based on quantiles of the data and is the
same in patsy and mgcv. Boundary points are at the limits of the data
range.
The available options used with `get_knots_bsplines` are
- knots : None or array
interior knots
- spacing : 'quantile' or 'equal'
- lower_bound : None or float
location of lower boundary knots, all boundary knots are at the same
point
- upper_bound : None or float
location of upper boundary knots, all boundary knots are at the same
point
- all_knots : None or array
If all knots are provided, then those will be taken as given and
all other options will be ignored.
"""
def __init__(self, x, df, degree=3, include_intercept=False,
constraints=None, variable_name="x",
covder2_kwds=None, **knot_kwds):
self.degree = degree
self.df = df
self.include_intercept = include_intercept
self.knots = get_knots_bsplines(x, degree=degree, df=df, **knot_kwds)
self.covder2_kwds = (covder2_kwds if covder2_kwds is not None
else {})
super().__init__(
x, constraints=constraints, variable_name=variable_name)
def _smooth_basis_for_single_variable(self):
basis, der_basis, der2_basis = _eval_bspline_basis(
self.x, self.knots, self.degree,
include_intercept=self.include_intercept)
# cov_der2 = np.dot(der2_basis.T, der2_basis)
cov_der2 = get_covder2(self, skip_ctransf=True,
**self.covder2_kwds)
return basis, der_basis, der2_basis, cov_der2
def transform(self, x_new, deriv=0, skip_ctransf=False):
"""
Create the spline basis for new observations
The main use of this stateful transformation is for prediction
using the same specification of the spline basis.
Parameters
----------
x_new : ndarray
observations of the underlying explanatory variable
deriv : int
which derivative of the spline basis to compute
This is an option for internal computation.
skip_ctransf : bool
whether to skip the constraint transform
This is an option for internal computation.
Returns
-------
basis : ndarray
design matrix for the spline basis for given ``x_new``
"""
if x_new is None:
x_new = self.x
exog = _eval_bspline_basis(x_new, self.knots, self.degree,
deriv=deriv,
include_intercept=self.include_intercept)
# ctransf does not exist yet when cov_der2 is computed
ctransf = getattr(self, "ctransf", None)
if ctransf is not None and not skip_ctransf:
exog = exog.dot(self.ctransf)
return exog
class UnivariateCubicSplines(UnivariateGamSmoother):
"""
Cubic Spline single smooth component
Cubic splines as described in Wood's book in chapter 3.
Parameters
----------
x : ndarray, 1-D
Underlying explanatory variable for the smooth term.
df : int
Number of basis functions or degrees of freedom.
constraints : {None, str, array_like}
Constraints are used to transform the basis functions to satisfy
those constraints. ``constraints = 'no-const'`` drops the
constant column from the basis.
transform : {None, "domain", tuple}
Rescaling applied to `x` before computing the basis. If
``"domain"``, `x` is rescaled to the unit interval using its own
minimum and maximum. If a tuple ``(low, upp)``, `x` is rescaled
to the unit interval using those bounds instead.
variable_name : str
The name for the underlying explanatory variable, x, used for
creating the column and parameter names for the basis functions.
"""
def __init__(self, x, df, constraints=None, transform="domain",
variable_name="x"):
self.degree = 3
self.df = df
self.transform_data_method = transform
self.x = x = self.transform_data(x, initialize=True)
self.knots = _equally_spaced_knots(x, df)
super().__init__(
x, constraints=constraints, variable_name=variable_name)
def transform_data(self, x, initialize=False):
"""
Rescale x according to `transform_data_method`
Parameters
----------
x : ndarray
Observations of the underlying explanatory variable.
initialize : bool
If True, the rescaling bounds (`domain_low`, `domain_upp`,
`domain_diff`) are computed from `x` (or from the tuple
passed to `transform` at instance construction) and stored
for later use.
Returns
-------
x : ndarray
Rescaled `x`.
"""
tm = self.transform_data_method
if tm is None:
return x
if initialize is True:
if tm == "domain":
self.domain_low = x.min(0)
self.domain_upp = x.max(0)
elif isinstance(tm, tuple):
self.domain_low = tm[0]
self.domain_upp = tm[1]
self.transform_data_method = "domain"
else:
raise ValueError("transform should be None, 'domain' "
"or a tuple")
self.domain_diff = self.domain_upp - self.domain_low
if self.transform_data_method == "domain":
x = (x - self.domain_low) / self.domain_diff
return x
else:
raise ValueError("incorrect transform_data_method")
def _smooth_basis_for_single_variable(self):
basis = self._splines_x()[:, :-1]
# demean except for constant, does not affect derivatives
if not self.constraints == "none":
self.transf_mean = basis[:, 1:].mean(0)
basis[:, 1:] -= self.transf_mean
else:
self.transf_mean = np.zeros(basis.shape[1])
s = self._splines_s()[:-1, :-1]
if not self.constraints == "none":
ctransf = np.diag(1/np.max(np.abs(basis), axis=0))
else:
ctransf = np.eye(basis.shape[1])
# use np.eye to avoid rescaling
# ctransf = np.eye(basis.shape[1])
if self.constraints == "no-const":
ctransf = ctransf[1:]
self.ctransf = ctransf
return basis, None, None, s
def _rk(self, x, z):
"""Cubic-spline reproducing-kernel function evaluated at x and z."""
p1 = ((z - 1 / 2) ** 2 - 1 / 12) * ((x - 1 / 2) ** 2 - 1 / 12) / 4
p2 = ((np.abs(z - x) - 1 / 2) ** 4 -
1 / 2 * (np.abs(z - x) - 1 / 2) ** 2 +
7 / 240) / 24.
return p1 - p2
def _splines_x(self, x=None):
"""
Build the cubic-spline design matrix at x, using self.knots
Parameters
----------
x : ndarray, optional
Observations of the underlying explanatory variable. If not
provided, uses `self.x`.
Returns
-------
basis : ndarray
Design matrix with an intercept column, a linear column, and
one column of `_rk` values for each knot.
"""
if x is None:
x = self.x
n_columns = len(self.knots) + 2
nobs = x.shape[0]
basis = np.ones(shape=(nobs, n_columns))
basis[:, 1] = x
# for loop equivalent to outer(x, xk, fun=rk)
for i, xi in enumerate(x):
for j, xkj in enumerate(self.knots):
s_ij = self._rk(xi, xkj)
basis[i, j + 2] = s_ij
return basis
def _splines_s(self):
"""Penalty matrix based on `_rk` evaluated at pairs of knots."""
q = len(self.knots) + 2
s = np.zeros(shape=(q, q))
for i, x1 in enumerate(self.knots):
for j, x2 in enumerate(self.knots):
s[i + 2, j + 2] = self._rk(x1, x2)
return s
def transform(self, x_new):
"""
Create the spline basis for new observations
Parameters
----------
x_new : ndarray
Observations of the underlying explanatory variable.
Returns
-------
exog : ndarray
Design matrix for the spline basis for given `x_new`.
"""
x_new = self.transform_data(x_new, initialize=False)
exog = self._splines_x(x_new)
exog[:, 1:] -= self.transf_mean
if self.ctransf is not None:
exog = exog.dot(self.ctransf)
return exog
class UnivariateCubicCyclicSplines(UnivariateGamSmoother):
"""
cyclic cubic regression spline single smooth component
This creates and holds the Cyclic CubicSpline basis function for one
component.
Parameters
----------
x : ndarray, 1-D
underlying explanatory variable for smooth terms.
df : int
number of basis functions or degrees of freedom
constraints : {None, str, array}
Constraints are used to transform the basis functions to satisfy
those constraints.
`constraints = 'center'` applies a linear transform to remove the
constant and center the basis functions.
variable_name : None or str
The name for the underlying explanatory variable, x, used in for
creating the column and parameter names for the basis functions.
"""
def __init__(self, x, df, constraints=None, variable_name="x"):
self.degree = 3
self.df = df
self.x = x
self.knots = _equally_spaced_knots(x, df)
super().__init__(
x, constraints=constraints, variable_name=variable_name)
def _smooth_basis_for_single_variable(self):
"""
Build the cyclic cubic regression spline basis and penalty
Returns
-------
basis : ndarray
Design matrix for the cyclic cubic spline basis.
der_basis : None
Not computed for this smoother.
der2_basis : None
Not computed for this smoother.
s : ndarray
Penalty matrix based on the second derivative.
"""
mgr = FormulaManager()
basis = mgr.get_matrices(
"cc(x, df=" + str(self.df) + ") - 1",
pd.DataFrame({"x": self.x}),
pandas=False,
)
self.model_spec = mgr.spec
n_inner_knots = self.df - 2 + 1 # +n_constraints
# TODO: from CubicRegressionSplines class
all_knots = get_all_sorted_knots(
self.x,
n_inner_knots=n_inner_knots,
inner_knots=None,
lower_bound=None,
upper_bound=None,
)
b, d = self._get_b_and_d(all_knots)
s = self._get_s(b, d)
return basis, None, None, s
def _get_b_and_d(self, knots):
"""
Returns mapping of cyclic cubic spline values to 2nd derivatives
.. note:: See 'Generalized Additive Models', Simon N. Wood, 2006,
pp 146-147
Parameters
----------
knots : ndarray
The 1-d array knots used for cubic spline parametrization,
must be sorted in ascending order.
Returns
-------
b : ndarray
Array for mapping cyclic cubic spline values at knots to
second derivatives.
d : ndarray
Array for mapping cyclic cubic spline values at knots to
second derivatives.
Notes
-----
The penalty matrix is equal to ``s = d.T.dot(b^-1).dot(d)``
"""
h = knots[1:] - knots[:-1]
n = knots.size - 1
# b and d are defined such that the penalty matrix is equivalent to:
# s = d.T.dot(b^-1).dot(d)
# reference in particular to pag 146 of Wood's book
b = np.zeros((n, n)) # the b matrix on page 146 of Wood's book
d = np.zeros((n, n)) # the d matrix on page 146 of Wood's book
b[0, 0] = (h[n - 1] + h[0]) / 3.
b[0, n - 1] = h[n - 1] / 6.
b[n - 1, 0] = h[n - 1] / 6.
d[0, 0] = -1. / h[0] - 1. / h[n - 1]
d[0, n - 1] = 1. / h[n - 1]
d[n - 1, 0] = 1. / h[n - 1]
for i in range(1, n):
b[i, i] = (h[i - 1] + h[i]) / 3.
b[i, i - 1] = h[i - 1] / 6.
b[i - 1, i] = h[i - 1] / 6.
d[i, i] = -1. / h[i - 1] - 1. / h[i]
d[i, i - 1] = 1. / h[i - 1]
d[i - 1, i] = 1. / h[i - 1]
return b, d
def _get_s(self, b, d):
"""Penalty matrix ``s = d.T.dot(b^-1).dot(d)``."""
return d.T.dot(np.linalg.inv(b)).dot(d)
def transform(self, x_new):
"""
Create the spline basis for new observations
Parameters
----------
x_new : ndarray
Observations of the underlying explanatory variable.
Returns
-------
exog : ndarray
Design matrix for the spline basis for given `x_new`.
"""
mgr = FormulaManager()
exog = mgr.get_matrices(
self.model_spec, pd.DataFrame({"x": x_new}), pandas=False
)
if self.ctransf is not None:
exog = exog.dot(self.ctransf)
return exog
class AdditiveGamSmoother(with_metaclass(ABCMeta)):
"""
Base class for additive smooth components
Parameters
----------
x : array_like, 1-D or 2-D
Underlying explanatory variable(s) for the smooth terms. If
2-dimensional, observations are in rows and explanatory
variables in columns.
variable_names : {list[str], None}
Names for the underlying explanatory variables. If None and
`x` is a pandas object, the names are taken from it, otherwise
generic names are generated.
include_intercept : bool or list[bool]
If False, then the basis functions of each smooth term are
transformed so that they do not include a constant. Can be a
single bool applied to all smooth terms, or a list with one
value per smooth term.
**kwargs
Additional keyword arguments passed through to subclasses.
"""
def __init__(self, x, variable_names=None, include_intercept=False,
**kwargs):
# get pandas names before using asarray
if isinstance(x, pd.DataFrame):
data_names = x.columns.tolist()
elif isinstance(x, pd.Series):
data_names = [x.name]
else:
data_names = None
x = np.asarray(x)
if x.ndim == 1:
self.x = x.copy()
self.x = inplace_reshape(self.x, (len(x), 1))
else:
self.x = x
self.nobs, self.k_variables = self.x.shape
if isinstance(include_intercept, bool):
self.include_intercept = [include_intercept] * self.k_variables
else:
self.include_intercept = include_intercept
if variable_names is None:
if data_names is not None:
self.variable_names = data_names
else:
self.variable_names = ["x" + str(i)
for i in range(self.k_variables)]
else:
self.variable_names = variable_names
self.smoothers = self._make_smoothers_list()
self.basis = np.hstack([smoother.basis for smoother in self.smoothers])
self.dim_basis = self.basis.shape[1]
self.penalty_matrices = [smoother.cov_der2
for smoother in self.smoothers]
self.col_names = []
for smoother in self.smoothers:
self.col_names.extend(smoother.col_names)
self.mask = []
last_column = 0
for smoother in self.smoothers:
mask = np.array([False] * self.dim_basis)
mask[last_column:smoother.dim_basis + last_column] = True
last_column = last_column + smoother.dim_basis
self.mask.append(mask)
@abstractmethod
def _make_smoothers_list(self):
pass
def transform(self, x_new):
"""
Create the spline basis for new observations
The main use of this stateful transformation is for prediction
using the same specification of the spline basis.
Parameters
----------
x_new: ndarray
observations of the underlying explanatory variable
Returns
-------
basis : ndarray
design matrix for the spline basis for given ``x_new``.
"""
if x_new.ndim == 1 and self.k_variables == 1:
x_new = x_new.reshape(-1, 1)
exog = np.hstack(
[self.smoothers[i].transform(x_new[:, i])
for i in range(self.k_variables)]
)
return exog
class GenericSmoothers(AdditiveGamSmoother):
"""
Generic class for additive smooth components for GAM
Parameters
----------
x : array_like, 1-D or 2-D
Underlying explanatory variable(s) for the smooth terms. If
2-dimensional, observations are in rows and explanatory
variables in columns.
smoothers : list of univariate smooth component instances
The individual smooth components making up the additive model.
"""
def __init__(self, x, smoothers):
self.smoothers = smoothers
super().__init__(x, variable_names=None)
def _make_smoothers_list(self):
return self.smoothers
class PolynomialSmoother(AdditiveGamSmoother):
"""
Additive polynomial components for GAM
Parameters
----------
x : array_like, 1-D or 2-D
Underlying explanatory variable(s) for the smooth terms. If
2-dimensional, observations are in rows and explanatory
variables in columns.
degrees : array_like[int]
Highest power included in the polynomial basis for each
explanatory variable in `x`.
variable_names : {list[str], None}
The names for the underlying explanatory variables, used for
creating the column and parameter names for the basis functions.
If ``x`` is a pandas object, then the names will be taken from it.
"""
def __init__(self, x, degrees, variable_names=None):
self.degrees = degrees
super().__init__(x, variable_names=variable_names)
def _make_smoothers_list(self):
smoothers = []
for v in range(self.k_variables):
uv_smoother = UnivariatePolynomialSmoother(
self.x[:, v],
degree=self.degrees[v],
variable_name=self.variable_names[v])
smoothers.append(uv_smoother)
return smoothers
[docs]
class BSplines(AdditiveGamSmoother):
"""
additive smooth components using B-Splines
This creates and holds the B-Spline basis function for several
components.
Parameters
----------
x : array_like, 1-D or 2-D
underlying explanatory variable for smooth terms.
If 2-dimensional, then observations should be in rows and
explanatory variables in columns.
df : {int, array_like[int]}
number of basis functions or degrees of freedom; should be equal
in length to the number of columns of `x`; may be an integer if
`x` has one column or is 1-D.
degree : {int, array_like[int]}
degree(s) of the spline; the same length and type rules apply as
to `df`
include_intercept : bool
If False, then the basis functions are transformed so that they
do not include a constant. This avoids perfect collinearity if
a constant or several components are included in the model.
constraints : {None, str, array}
Constraints are used to transform the basis functions to satisfy
those constraints.
`constraints = 'center'` applies a linear transform to remove the
constant and center the basis functions.
variable_names : {list[str], None}
The names for the underlying explanatory variables, x used in for
creating the column and parameter names for the basis functions.
If ``x`` is a pandas object, then the names will be taken from it.
knot_kwds : None or list of dict
option for the knot selection.
By default knots are selected in the same way as in patsy, however the
number of knots is independent of keeping or removing the constant.
Interior knot selection is based on quantiles of the data and is the
same in patsy and mgcv. Boundary points are at the limits of the data
range.
The available options used with `get_knots_bsplines` are
- knots : None or array
interior knots
- spacing : 'quantile' or 'equal'
- lower_bound : None or float
location of lower boundary knots, all boundary knots are at the same
point
- upper_bound : None or float
location of upper boundary knots, all boundary knots are at the same
point
- all_knots : None or array
If all knots are provided, then those will be taken as given and
all other options will be ignored.
Attributes
----------
smoothers : list of univariate smooth component instances
basis : design matrix, array of spline bases columns for all components
penalty_matrices : list of penalty matrices, one for each smooth term
dim_basis : number of columns in the basis
k_variables : number of smooth components
col_names : created names for the basis columns
There are additional attributes about the specification of the splines
and some attributes mainly for internal use.
Notes
-----
A constant in the spline basis function can be removed in two different
ways.
The first is by dropping one basis column and normalizing the
remaining columns. This is obtained by the default
``include_intercept=False, constraints=None``
The second option is by using the centering transform which is a linear
transformation of all basis functions. As a consequence of the
transformation, the B-spline basis functions do not have locally bounded
support anymore. This is obtained ``constraints='center'``. In this case
``include_intercept`` will be automatically set to True to avoid
dropping an additional column.
"""
def __init__(self, x, df, degree, include_intercept=False,
constraints=None, variable_names=None, knot_kwds=None):
if isinstance(degree, int):
self.degrees = np.array([degree], dtype=int)
else:
self.degrees = degree
if isinstance(df, int):
self.dfs = np.array([df], dtype=int)
else:
self.dfs = df
self.knot_kwds = knot_kwds
# TODO: move attaching constraints to super call
self.constraints = constraints
if constraints == "center":
include_intercept = True
super().__init__(x, include_intercept=include_intercept,
variable_names=variable_names)
def _make_smoothers_list(self):
smoothers = []
for v in range(self.k_variables):
kwds = self.knot_kwds[v] if self.knot_kwds else {}
uv_smoother = UnivariateBSplines(
self.x[:, v],
df=self.dfs[v], degree=self.degrees[v],
include_intercept=self.include_intercept[v],
constraints=self.constraints,
variable_name=self.variable_names[v], **kwds)
smoothers.append(uv_smoother)
return smoothers
class CubicSplines(AdditiveGamSmoother):
"""
additive smooth components using cubic splines as in Wood 2006
Note, these splines do NOT use the same spline basis as
``Cubic Regression Splines``.
Parameters
----------
x : array_like, 1-D or 2-D
underlying explanatory variable for smooth terms.
If 2-dimensional, then observations should be in rows and
explanatory variables in columns.
df : array_like[int]
number of basis functions or degrees of freedom; should be equal
in length to the number of columns of `x`.
constraints : {None, str, array}
Constraints are used to transform the basis functions to satisfy
those constraints.
`constraints = 'center'` applies a linear transform to remove the
constant and center the basis functions.
transform : {None, "domain", tuple}
Rescaling applied to each column of `x` before computing the
basis. See `UnivariateCubicSplines` for details.
variable_names : {list[str], None}
The names for the underlying explanatory variables, x used in for
creating the column and parameter names for the basis functions.
If ``x`` is a pandas object, then the names will be taken from it.
"""
def __init__(self, x, df, constraints="center", transform="domain",
variable_names=None):
self.dfs = df
self.constraints = constraints
self.transform = transform
super().__init__(x, constraints=constraints,
variable_names=variable_names)
def _make_smoothers_list(self):
smoothers = []
for v in range(self.k_variables):
uv_smoother = UnivariateCubicSplines(
self.x[:, v], df=self.dfs[v],
constraints=self.constraints,
transform=self.transform,
variable_name=self.variable_names[v])
smoothers.append(uv_smoother)
return smoothers
[docs]
class CyclicCubicSplines(AdditiveGamSmoother):
"""
additive smooth components using cyclic cubic regression splines
This spline basis is the same as in patsy.
Parameters
----------
x : array_like, 1-D or 2-D
underlying explanatory variable for smooth terms.
If 2-dimensional, then observations should be in rows and
explanatory variables in columns.
df : int
number of basis functions or degrees of freedom
constraints : {None, str, array}
Constraints are used to transform the basis functions to satisfy
those constraints.
variable_names : {list[str], None}
The names for the underlying explanatory variables, x used in for
creating the column and parameter names for the basis functions.
If ``x`` is a pandas object, then the names will be taken from it.
"""
def __init__(self, x, df, constraints=None, variable_names=None):
self.dfs = df
# TODO: move attaching constraints to super call
self.constraints = constraints
super().__init__(x, variable_names=variable_names)
def _make_smoothers_list(self):
smoothers = []
for v in range(self.k_variables):
uv_smoother = UnivariateCubicCyclicSplines(
self.x[:, v],
df=self.dfs[v], constraints=self.constraints,
variable_name=self.variable_names[v])
smoothers.append(uv_smoother)
return smoothers
# class CubicRegressionSplines(BaseCubicSplines):
# # TODO: this class is still not tested
#
# def __init__(self, x, df=10):
# import warnings
# warnings.warn("This class is still not tested and it is probably"
# " not working properly. "
# "I suggest to use another smoother", Warning)
#
# super(CubicRegressionSplines, self).__init__(x, df)
#
# self.basis = dmatrix("cc(x, df=" + str(df) + ") - 1", {"x": x})
# n_inner_knots = df - 2 + 1 # +n_constraints
# # TODO: ACcording to CubicRegressionSplines class this should be
# # n_inner_knots = df - 2
# all_knots = _get_all_sorted_knots(x, n_inner_knots=n_inner_knots,
# inner_knots=None,
# lower_bound=None, upper_bound=None)
#
# b, d = self._get_b_and_d(all_knots)
# self.s = self._get_s(b, d)
#
# self.dim_basis = self.basis.shape[1]
#
# def _get_b_and_d(self, knots):
#
# h = knots[1:] - knots[:-1]
# n = knots.size - 1
#
# # b and d are defined such that the penalty matrix is equivalent to:
# # s = d.T.dot(b^-1).dot(d)
# # reference in particular to pag 146 of Wood's book
# b = np.zeros((n, n)) # the b matrix on page 146 of Wood's book
# d = np.zeros((n, n)) # the d matrix on page 146 of Wood's book
#
# for i in range(n-2):
# d[i, i] = 1/h[i]
# d[i, i+1] = -1/h[i] - 1/h[i+1]
# d[i, i+2] = 1/h[i+1]
#
# b[i, i] = (h[i] + h[i+1])/3
#
# for i in range(n-3):
# b[i, i+1] = h[i+1]/6
# b[i+1, i] = h[i+1]/6
#
# return b, d
#
# def _get_s(self, b, d):
#
# return d.T.dot(np.linalg.pinv(b)).dot(d)
