Welcome to py-earth’s documentation!

Introduction

The py-earth package is a Python implementation of Jerome Friedman’s Multivariate Adaptive Regression Splines algorithm, in the style of scikit-learn. For more information about Multivariate Adaptive Regression Splines, see below. Py-earth is written in Python and Cython. It provides an interface that is compatible with scikit-learn’s Estimator, Predictor, Transformer, and Model interfaces. Py-earth accommodates input in the form of numpy arrays, pandas DataFrames, patsy DesignMatrix objects, or most anything that can be converted into an arrray of floats. Fitted models can be pickled for later use.

Multivariate Adaptive Regression Splines

Multivariate adaptive regression splines, implemented by the Earth class, is a flexible regression method that automatically searches for interactions and non-linear relationships. Earth models can be thought of as linear models in a higher dimensional basis space. Each term in an Earth model is a product of so called “hinge functions”. A hinge function is a function that’s equal to its argument where that argument is greater than zero and is zero everywhere else.

\[\begin{split}\text{h}\left(x-t\right)=\left[x-t\right]_{+}=\begin{cases} x-t, & x>t\\ 0, & x\leq t \end{cases}\end{split}\]
_images/hinge.png

An Earth model is a linear combination of basis functions, each of which is a product of one or more of the following:

  1. A constant
  2. Linear functions of input variables
  3. Hinge functions of input variables

For example, a simple piecewise linear function in one variable can be expressed as a linear combination of two hinge functions and a constant (see below). During fitting, the Earth class automatically determines which variables and basis functions to use. The algorithm has two stages. First, the forward pass searches for terms that locally minimize squared error loss on the training set. Next, a pruning pass selects a subset of those terms that produces a locally minimal generalized cross-validation (GCV) score. The GCV score is not actually based on cross-validation, but rather is meant to approximate a true cross-validation score by penalizing model complexity. The final result is a set of basis functions that is nonlinear in the original feature space, may include interactions, and is likely to generalize well.

\[y=1-2\text{h}\left(1-x\right)+\frac{1}{2}\text{h}\left(x-1\right)\]
_images/piecewise_linear.png

A Simple Earth Example

import numpy
from pyearth import Earth
from matplotlib import pyplot

#Create some fake data
numpy.random.seed(0)
m = 1000
n = 10
X = 80*numpy.random.uniform(size=(m,n)) - 40
y = numpy.abs(X[:,6] - 4.0) + 1*numpy.random.normal(size=m)

#Fit an Earth model
model = Earth()
model.fit(X,y)

#Print the model
print model.trace()
print model.summary()

#Plot the model
y_hat = model.predict(X)
pyplot.figure()
pyplot.plot(X[:,6],y,'r.')
pyplot.plot(X[:,6],y_hat,'b.')
pyplot.xlabel('x_6')
pyplot.ylabel('y')
pyplot.title('Simple Earth Example')
pyplot.show()
_images/simple_earth_example.png

Bibliography

[Bjo96]Ake Bjorck. Numerical Methods for Least Squares Problems. Society for Industrial and Applied Mathematics, Philadelphia, 1996.
[Fri93](1, 2) Jerome H. Friedman. Technical Report No. 110: Fast MARS.. Technical Report, Stanford University Department of Statistics, 1993. URL: http://scholar.google.com/scholar?hl=en&btnG=Search&q=intitle:Fast+MARS#0.
[Fri91a](1, 2) JH Friedman. Multivariate adaptive regression splines. The annals of statistics, 19(1):1–67, 1991. URL: http://www.jstor.org/stable/10.2307/2241837.
[Fri91b](1, 2) JH Friedman. Technical Report No. 108: Estimating functions of mixed ordinal and categorical variables using adaptive splines. Technical Report, Stanford University Department of Statistics, 1991. URL: http://scholar.google.com/scholar?hl=en&btnG=Search&q=intitle:Estimating+functions+of+mixed+ordinal+and+categorical+variables+using+adaptive+splines#0.
[GVanLoan96]Gene Golub and Charles Van Loan. Matrix Computations. Johns Hopkins University Press, 3 edition, 1996.
[HTF09]Trevor Hastie, Robert Tibshirani, and Jerome Friedman. Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer Science+Business Media, New York, 2 edition, 2009.
[Mil12](1, 2) Stephen Millborrow. earth: Multivariate Adaptive Regression Spline Models. 2012. URL: http://cran.r-project.org/web/packages/earth/index.html.
[Ste98]G. W. Stewart. Matrix Algorithms Volume 1: Basic Decompositions. Society for Industrial and Applied Mathematics, Philadelphia, 1998.

References [HTF09], [Mil12], [Fri91a], [Fri93], and [Fri91b] contain discussions likely to be useful to users of py-earth. References [Fri91a], [Mil12], [Bjo96], [Ste98], [GVanLoan96], [Fri93], and [Fri91b] were useful during the implementation process.

API

class pyearth.Earth(max_terms=None, max_degree=None, penalty=None, endspan_alpha=None, endspan=None, minspan_alpha=None, minspan=None, thresh=None, min_search_points=None, check_every=None, allow_linear=None, smooth=None, enable_pruning=True)[source]

Multivariate Adaptive Regression Splines

A flexible regression method that automatically searches for interactions and non-linear relationships. Earth models can be thought of as linear models in a higher dimensional basis space (specifically, a multivariate truncated power spline basis). Each term in an Earth model is a product of so called “hinge functions”. A hinge function is a function that’s equal to its argument where that argument is greater than zero and is zero everywhere else.

The multivariate adaptive regression splines algorithm has two stages. First, the forward pass searches for terms in the truncated power spline basis that locally minimize the squared error loss of the training set. Next, a pruning pass selects a subset of those terms that produces a locally minimal generalized cross-validation (GCV) score. The GCV score is not actually based on cross-validation, but rather is meant to approximate a true cross-validation score by penalizing model complexity. The final result is a set of terms that is nonlinear in the original feature space, may include interactions, and is likely to generalize well.

The Earth class supports dense input only. Data structures from the pandas and patsy modules are supported, but are copied into numpy arrays for computation. No copy is made if the inputs are numpy float64 arrays. Earth objects can be serialized using the pickle module and copied using the copy module.

Parameters:

max_terms : int, optional (default=2*n + 10, where n is the

number of features)

The maximum number of terms generated by the forward pass.

max_degree : int, optional (default=1)

The maximum degree of terms generated by the forward pass.

penalty : float, optional (default=3.0)

A smoothing parameter used to calculate GCV and GRSQ. Used during the pruning pass and to determine whether to add a hinge or linear basis function during the forward pass. See the d parameter in equation 32, Friedman, 1991.

endspan_alpha : float, optional, probability between 0 and 1 (default=0.05)

A parameter controlling the calculation of the endspan parameter (below). The endspan parameter is calculated as round(3 - log2(endspan_alpha/n)), where n is the number of features. The endspan_alpha parameter represents the probability of a run of positive or negative error values on either end of the data vector of any feature in the data set. See equation 45, Friedman, 1991.

endspan : int, optional (default=-1)

The number of extreme data values of each feature not eligible as knot locations. If endspan is set to -1 (default) then the endspan parameter is calculated based on endspan_alpah (above). If endspan is set to a positive integer then endspan_alpha is ignored.

minspan_alpha : float, optional, probability between 0 and 1 (default=0.05)

A parameter controlling the calculation of the minspan parameter (below). The minspan parameter is calculated as

(int) -log2(-(1.0/(n*count))*log(1.0-minspan_alpha)) / 2.5

where n is the number of features and count is the number of points at which the parent term is non-zero. The minspan_alpha parameter represents the probability of a run of positive or negative error values between adjacent knots separated by minspan intervening data points. See equation 43, Friedman, 1991.

minspan : int, optional (default=-1)

The minimal number of data points between knots. If minspan is set to -1 (default) then the minspan parameter is calculated based on minspan_alpha (above). If minspan is set to a positive integer then minspan_alpha is ignored.

thresh : float, optional (defaul=0.001)

Parameter used when evaluating stopping conditions for the forward pass. If either RSQ > 1 - thresh or if RSQ increases by less than thresh for a forward pass iteration then the forward pass is terminated.

min_search_points : int, optional (default=100)

Used to calculate check_every (below). The minimum samples necessary for check_every to be greater than 1. The check_every parameter is calculated as

(int) m / min_search_points

if m > min_search_points, where m is the number of samples in the training set. If m <= min_search_points then check_every is set to 1.

check_every : int, optional (default=-1)

If check_every > 0, only one of every check_every sorted data points is considered as a candidate knot. If check_every is set to -1 then the check_every parameter is calculated based on min_search_points (above).

allow_linear : bool, optional (default=True)

If True, the forward pass will check the GCV of each new pair of terms and, if it’s not an improvement on a single term with no knot (called a linear term, although it may actually be a product of a linear term with some other parent term), then only that single, knotless term will be used. If False, that behavior is disabled and all terms will have knots except those with variables specified by the linvars argument (see the fit method).

smooth : bool, optional (default=False)

If True, the model will be smoothed such that it has continuous first derivatives. For details, see section 3.7, Friedman, 1991.

enable_pruning : bool, optional(default=True)

If False, the pruning pass will be skipped.

References

[R1]Friedman, Jerome. Multivariate Adaptive Regression Splines. Annals of Statistics. Volume 19, Number 1 (1991), 1-67.

Attributes

coef_ (array, shape = [pruned basis length]) The weights of the model terms that have not been pruned.
basis_ (_basis.Basis) An object representing model terms. Each term is a product of constant, linear, and hinge functions of the input features.
mse_ (float) The mean squared error of the model after the final linear fit. If sample_weight is given, this score is weighted appropriately.
rsq_ (float) The generalized r^2 of the model after the final linear fit. If sample_weight is given, this score is weighted appropriately.
gcv_ (float) The generalized cross validation (GCV) score of the model after the final linear fit. If sample_weight is given, this score is weighted appropriately.
grsq_ (float) An r^2 like score based on the GCV. If sample_weight is given, this score is weighted appropriately.
forward_pass_record_ (_record.ForwardPassRecord) An object containing information about the forward pass, such as training loss function values after each iteration and the final stopping condition.
pruning_pass_record_ (_record.PruningPassRecord) An object containing information about the pruning pass, such as training loss function values after each iteration and the selected optimal iteration.
xlabels_ (list) List of column names for training predictors. Defaults to [‘x0’,’x1’,....] if column names are not provided.
fit(X, y=None, sample_weight=None, xlabels=None, linvars=[])[source]

Fit an Earth model to the input data X and y.

Parameters:

X : array-like, shape = [m, n] where m is the number of samples

and n is the number of features The training predictors. The X parameter can be a numpy array, a pandas DataFrame, a patsy DesignMatrix, or a tuple of patsy DesignMatrix objects as output by patsy.dmatrices.

y : array-like, optional (default=None), shape = [m] where m is the

number of samples The training response. The y parameter can be a numpy array, a pandas DataFrame with one column, a Patsy DesignMatrix, or can be left as None (default) if X was the output of a call to patsy.dmatrices (in which case, X contains the response).

sample_weight : array-like, optional (default=None), shape = [m] where

m is the number of samples Sample weights for training. Weights must be greater than or equal to zero. Rows with greater weights contribute more strongly to the fitted model. Rows with zero weight do not contribute at all. Weights are useful when dealing with heteroscedasticity. In such cases, the weight should be proportional to the inverse of the (known) variance.

linvars : iterable of strings or ints, optional (empty by default)

Used to specify features that may only enter terms as linear basis functions (without knots). Can include both column numbers and column names (see xlabels, below). If left empty, some variables may still enter linearly during the forward pass if no knot would provide a reduction in GCV compared to the linear function. Note that this feature differs from the R package earth.

xlabels : iterable of strings, optional (empty by default)

The xlabels argument can be used to assign names to data columns. This argument is not generally needed, as names can be captured automatically from most standard data structures. If included, must have length n, where n is the number of features. Note that column order is used to compute term values and make predictions, not column names.

fit_transform(X, y=None, **fit_params)

Fit to data, then transform it.

Fits transformer to X and y with optional parameters fit_params and returns a transformed version of X.

Parameters:

X : numpy array of shape [n_samples, n_features]

Training set.

y : numpy array of shape [n_samples]

Target values.

Returns:

X_new : numpy array of shape [n_samples, n_features_new]

Transformed array.

forward_trace()[source]

Return information about the forward pass.

get_params(deep=True)

Get parameters for this estimator.

Parameters:

deep: boolean, optional

If True, will return the parameters for this estimator and contained subobjects that are estimators.

Returns:

params : mapping of string to any

Parameter names mapped to their values.

get_penalty()[source]

Get the penalty parameter being used. Default is 3.

predict(X)[source]

Predict the response based on the input data X.

Parameters:

X : array-like, shape = [m, n] where m is the number of samples and n

is the number of features The training predictors. The X parameter can be a numpy array, a pandas DataFrame, or a patsy DesignMatrix.

predict_deriv(X, variables=None)[source]

Predict the first derivatives of the response based on the input data X.

Parameters:

X : array-like, shape = [m, n] where m is the number of samples and n is

the number of features The training predictors. The X parameter can be a numpy array, a pandas DataFrame, or a patsy DesignMatrix.

variables : list

The variables over which derivatives will be computed. Each column in the resulting array corresponds to a variable. If not specified, all variables are used (even if some are not relevant to the final model and have derivatives that are identically zero).

pruning_trace()[source]

Return information about the pruning pass.

score(X, y=None, sample_weight=None)[source]

Calculate the generalized r^2 of the model on data X and y.

Parameters:

X : array-like, shape = [m, n] where m is the number of samples

and n is the number of features The training predictors. The X parameter can be a numpy array, a pandas DataFrame, a patsy DesignMatrix, or a tuple of patsy DesignMatrix objects as output by patsy.dmatrices.

y : array-like, optional (default=None), shape = [m] where m is the

number of samples The training response. The y parameter can be a numpy array, a pandas DataFrame with one column, a Patsy DesignMatrix, or can be left as None (default) if X was the output of a call to patsy.dmatrices (in which case, X contains the response).

sample_weight : array-like, optional (default=None), shape = [m] where

m is the number of samples

Sample weights for training. Weights must be greater than or equal to zero. Rows with greater weights contribute more strongly to the fitted model. Rows with zero weight do not contribute at all. Weights are useful when dealing with heteroscedasticity. In such cases, the weight should be proportional to the inverse of the (known) variance.

set_params(**params)

Set the parameters of this estimator.

The method works on simple estimators as well as on nested objects (such as pipelines). The former have parameters of the form <component>__<parameter> so that it’s possible to update each component of a nested object.

Returns:self
summary()[source]

Return a string describing the model.

trace()[source]

Return information about the forward and pruning passes.

transform(X)[source]

Transform X into the basis space. Normally, users will call the predict method instead, which both transforms into basis space calculates the weighted sum of basis terms to produce a prediction of the response. Users may wish to call transform directly in some cases. For example, users may wish to apply other statistical or machine learning algorithms, such as generalized linear regression, in basis space.

Parameters:

X : array-like, shape = [m, n] where m is the number of samples and n

is the number of features The training predictors. The X parameter can be a numpy array, a pandas DataFrame, or a patsy DesignMatrix.