1.1.1.Ordinary Least Squares#

LinearRegression fits a linear model with coefficients\(w = (w_1, ..., w_p)\) to minimize the residual sumof squares between the observed targets in the dataset, and thetargets predicted by the linear approximation. Mathematically itsolves a problem of the form:

LinearRegression takes in itsfit method argumentsX,y,sample_weight and stores the coefficients\(w\) of the linear model in itscoef_ andintercept_ attributes:

>>>fromsklearnimportlinear_model>>>reg=linear_model.LinearRegression()>>>reg.fit([[0,0],[1,1],[2,2]],[0,1,2])LinearRegression()>>>reg.coef_array([0.5, 0.5])>>>reg.intercept_0.0

The coefficient estimates for Ordinary Least Squares rely on theindependence of the features. When features are correlated and somecolumns of the design matrix\(X\) have an approximately lineardependence, the design matrix becomes close to singularand as a result, the least-squares estimate becomes highly sensitiveto random errors in the observed target, producing a largevariance. This situation ofmulticollinearity can arise, forexample, when data are collected without an experimental design.

Examples

• Ordinary Least Squares and Ridge Regression

1.1.2.Ridge regression and classification#

1.1.2.1.Regression#

Ridge regression addresses some of the problems ofOrdinary Least Squares by imposing a penalty on the size of thecoefficients. The ridge coefficients minimize a penalized residual sumof squares:

\[\min_{w} || X w - y||_2^2 + \alpha ||w||_2^2\]

The complexity parameter\(\alpha \geq 0\) controls the amountof shrinkage: the larger the value of\(\alpha\), the greater the amountof shrinkage and thus the coefficients become more robust to collinearity.

As with other linear models,Ridge will take in itsfit methodarraysX,y and will store the coefficients\(w\) of the linear model initscoef_ member:

>>>fromsklearnimportlinear_model>>>reg=linear_model.Ridge(alpha=.5)>>>reg.fit([[0,0],[0,0],[1,1]],[0,.1,1])Ridge(alpha=0.5)>>>reg.coef_array([0.34545455, 0.34545455])>>>reg.intercept_np.float64(0.13636)

Note that the classRidge allows for the user to specify that thesolver be automatically chosen by settingsolver="auto". When this optionis specified,Ridge will choose between the"lbfgs","cholesky",and"sparse_cg" solvers.Ridge will begin checking the conditionsshown in the following table from top to bottom. If the condition is true,the corresponding solver is chosen.

Solver	Condition
‘lbfgs’	Thepositive=True option is specified.
‘cholesky’	The input array X is not sparse.
‘sparse_cg’	None of the above conditions are fulfilled.

Examples

• Ordinary Least Squares and Ridge Regression

• Plot Ridge coefficients as a function of the regularization

• Common pitfalls in the interpretation of coefficients of linear models

>>>importnumpyasnp>>>fromsklearnimportlinear_model>>>reg=linear_model.RidgeCV(alphas=np.logspace(-6,6,13))>>>reg.fit([[0,0],[0,0],[1,1]],[0,.1,1])RidgeCV(alphas=array([1.e-06, 1.e-05, 1.e-04, 1.e-03, 1.e-02, 1.e-01, 1.e+00, 1.e+01,      1.e+02, 1.e+03, 1.e+04, 1.e+05, 1.e+06]))>>>reg.alpha_np.float64(0.01)

1.1.3.Lasso#

TheLasso is a linear model that estimates sparse coefficients.It is useful in some contexts due to its tendency to prefer solutionswith fewer non-zero coefficients, effectively reducing the number offeatures upon which the given solution is dependent. For this reason,Lasso and its variants are fundamental to the field of compressed sensing.Under certain conditions, it can recover the exact set of non-zerocoefficients (seeCompressive sensing: tomography reconstruction with L1 prior (Lasso)).

Mathematically, it consists of a linear model with an added regularization term.The objective function to minimize is:

The lasso estimate thus solves the minimization of theleast-squares penalty with\(\alpha ||w||_1\) added, where\(\alpha\) is a constant and\(||w||_1\) is the\(\ell_1\)-norm ofthe coefficient vector.

The implementation in the classLasso uses coordinate descent asthe algorithm to fit the coefficients. SeeLeast Angle Regressionfor another implementation:

>>>fromsklearnimportlinear_model>>>reg=linear_model.Lasso(alpha=0.1)>>>reg.fit([[0,0],[1,1]],[0,1])Lasso(alpha=0.1)>>>reg.predict([[1,1]])array([0.8])

The functionlasso_path is useful for lower-level tasks, as itcomputes the coefficients along the full path of possible values.

Examples

• L1-based models for Sparse Signals

• Compressive sensing: tomography reconstruction with L1 prior (Lasso)

• Common pitfalls in the interpretation of coefficients of linear models

1.1.3.1.Setting regularization parameter#

Thealpha parameter controls the degree of sparsity of the estimatedcoefficients.

1.1.3.1.1.Using cross-validation#

scikit-learn exposes objects that set the Lassoalpha parameter bycross-validation:LassoCV andLassoLarsCV.LassoLarsCV is based on theLeast Angle Regression algorithmexplained below.

For high-dimensional datasets with many collinear features,LassoCV is most often preferable. However,LassoLarsCV hasthe advantage of exploring more relevant values ofalpha parameter, andif the number of samples is very small compared to the number offeatures, it is often faster thanLassoCV.

1.1.3.1.2.Information-criteria based model selection#

Alternatively, the estimatorLassoLarsIC proposes to use theAkaike information criterion (AIC) and the Bayes Information criterion (BIC).It is a computationally cheaper alternative to find the optimal value of alphaas the regularization path is computed only once instead of k+1 timeswhen using k-fold cross-validation.

Indeed, these criteria are computed on the in-sample training set. In short,they penalize the over-optimistic scores of the different Lasso models bytheir flexibility (cf. to “Mathematical details” section below).

However, such criteria need a proper estimation of the degrees of freedom ofthe solution, are derived for large samples (asymptotic results) and assume thecorrect model is candidates under investigation. They also tend to break whenthe problem is badly conditioned (e.g. more features than samples).

Examples

• Lasso model selection: AIC-BIC / cross-validation

• Lasso model selection via information criteria

1.1.3.1.3.AIC and BIC criteria#

The definition of AIC (and thus BIC) might differ in the literature. In thissection, we give more information regarding the criterion computed inscikit-learn.

Mathematical details#

The AIC criterion is defined as:

\[AIC = -2 \log(\hat{L}) + 2 d\]

where\(\hat{L}\) is the maximum likelihood of the model and\(d\) is the number of parameters (as well referred to as degrees offreedom in the previous section).

The definition of BIC replaces the constant\(2\) by\(\log(N)\):

\[BIC = -2 \log(\hat{L}) + \log(N) d\]

where\(N\) is the number of samples.

For a linear Gaussian model, the maximum log-likelihood is defined as:

\[\log(\hat{L}) = - \frac{n}{2} \log(2 \pi) - \frac{n}{2} \log(\sigma^2) - \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{2\sigma^2}\]

where\(\sigma^2\) is an estimate of the noise variance,\(y_i\) and\(\hat{y}_i\) are respectively the true and predictedtargets, and\(n\) is the number of samples.

Plugging the maximum log-likelihood in the AIC formula yields:

\[AIC = n \log(2 \pi \sigma^2) + \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{\sigma^2} + 2 d\]

The first term of the above expression is sometimes discarded since it is aconstant when\(\sigma^2\) is provided. In addition,it is sometimes stated that the AIC is equivalent to the\(C_p\) statistic[12]. In a strict sense, however, it is equivalent only up to some constantand a multiplicative factor.

At last, we mentioned above that\(\sigma^2\) is an estimate of thenoise variance. InLassoLarsIC when the parameternoise_variance isnot provided (default), the noise variance is estimated via the unbiasedestimator[13] defined as:

\[\sigma^2 = \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{n - p}\]

where\(p\) is the number of features and\(\hat{y}_i\) is thepredicted target using an ordinary least squares regression. Note, that thisformula is valid only whenn_samples>n_features.

References

[12]

Zou, Hui, Trevor Hastie, and Robert Tibshirani.“On the degrees of freedom of the lasso.”The Annals of Statistics 35.5 (2007): 2173-2192.

[13]

Cherkassky, Vladimir, and Yunqian Ma.“Comparison of model selection for regression.”Neural computation 15.7 (2003): 1691-1714.

1.1.3.1.4.Comparison with the regularization parameter of SVM#

The equivalence betweenalpha and the regularization parameter of SVM,C is given byalpha=1/C oralpha=1/(n_samples*C),depending on the estimator and the exact objective function optimized by themodel.

1.1.4.Multi-task Lasso#

TheMultiTaskLasso is a linear model that estimates sparsecoefficients for multiple regression problems jointly:y is a 2D array,of shape(n_samples,n_tasks). The constraint is that the selectedfeatures are the same for all the regression problems, also called tasks.

The following figure compares the location of the non-zero entries in thecoefficient matrix W obtained with a simple Lasso or a MultiTaskLasso.The Lasso estimates yield scattered non-zeros while the non-zeros ofthe MultiTaskLasso are full columns.

Fitting a time-series model, imposing that any active feature be active at all times.

Examples

• Joint feature selection with multi-task Lasso

1.1.5.Elastic-Net#

ElasticNet is a linear regression model trained with both\(\ell_1\) and\(\ell_2\)-norm regularization of the coefficients.This combination allows for learning a sparse model where few ofthe weights are non-zero likeLasso, while still maintainingthe regularization properties ofRidge. We control the convexcombination of\(\ell_1\) and\(\ell_2\) using thel1_ratioparameter.

Elastic-net is useful when there are multiple features that arecorrelated with one another. Lasso is likely to pick one of theseat random, while elastic-net is likely to pick both.

A practical advantage of trading-off between Lasso and Ridge is that itallows Elastic-Net to inherit some of Ridge’s stability under rotation.

The objective function to minimize is in this case

The classElasticNetCV can be used to set the parametersalpha (\(\alpha\)) andl1_ratio (\(\rho\)) by cross-validation.

Examples

• L1-based models for Sparse Signals

• Lasso, Lasso-LARS, and Elastic Net paths

• Fitting an Elastic Net with a precomputed Gram Matrix and Weighted Samples

1.1.6.Multi-task Elastic-Net#

TheMultiTaskElasticNet is an elastic-net model that estimates sparsecoefficients for multiple regression problems jointly:Y is a 2D arrayof shape(n_samples,n_tasks). The constraint is that the selectedfeatures are the same for all the regression problems, also called tasks.

Mathematically, it consists of a linear model trained with a mixed\(\ell_1\)\(\ell_2\)-norm and\(\ell_2\)-norm for regularization.The objective function to minimize is:

The implementation in the classMultiTaskElasticNet uses coordinate descent asthe algorithm to fit the coefficients.

The classMultiTaskElasticNetCV can be used to set the parametersalpha (\(\alpha\)) andl1_ratio (\(\rho\)) by cross-validation.

1.1.7.Least Angle Regression#

Least-angle regression (LARS) is a regression algorithm forhigh-dimensional data, developed by Bradley Efron, Trevor Hastie, IainJohnstone and Robert Tibshirani. LARS is similar to forward stepwiseregression. At each step, it finds the feature most correlated with thetarget. When there are multiple features having equal correlation, insteadof continuing along the same feature, it proceeds in a direction equiangularbetween the features.

The advantages of LARS are:

• It is numerically efficient in contexts where the number of featuresis significantly greater than the number of samples.

• It is computationally just as fast as forward selection and hasthe same order of complexity as ordinary least squares.

• It produces a full piecewise linear solution path, which isuseful in cross-validation or similar attempts to tune the model.

• If two features are almost equally correlated with the target,then their coefficients should increase at approximately the samerate. The algorithm thus behaves as intuition would expect, andalso is more stable.

• It is easily modified to produce solutions for other estimators,like the Lasso.

The disadvantages of the LARS method include:

• Because LARS is based upon an iterative refitting of theresiduals, it would appear to be especially sensitive to theeffects of noise. This problem is discussed in detail by Weisbergin the discussion section of the Efron et al. (2004) Annals ofStatistics article.

The LARS model can be used via the estimatorLars, or itslow-level implementationlars_path orlars_path_gram.

1.1.8.LARS Lasso#

LassoLars is a lasso model implemented using the LARSalgorithm, and unlike the implementation based on coordinate descent,this yields the exact solution, which is piecewise linear as afunction of the norm of its coefficients.

>>>fromsklearnimportlinear_model>>>reg=linear_model.LassoLars(alpha=.1)>>>reg.fit([[0,0],[1,1]],[0,1])LassoLars(alpha=0.1)>>>reg.coef_array([0.6, 0.        ])

Examples

• Lasso, Lasso-LARS, and Elastic Net paths

The LARS algorithm provides the full path of the coefficients alongthe regularization parameter almost for free, thus a common operationis to retrieve the path with one of the functionslars_pathorlars_path_gram.

1.1.9.Orthogonal Matching Pursuit (OMP)#

OrthogonalMatchingPursuit andorthogonal_mp implement the OMPalgorithm for approximating the fit of a linear model with constraints imposedon the number of non-zero coefficients (i.e. the\(\ell_0\) pseudo-norm).

Being a forward feature selection method likeLeast Angle Regression,orthogonal matching pursuit can approximate the optimum solution vector with afixed number of non-zero elements:

Alternatively, orthogonal matching pursuit can target a specific error insteadof a specific number of non-zero coefficients. This can be expressed as:

OMP is based on a greedy algorithm that includes at each step the atom mosthighly correlated with the current residual. It is similar to the simplermatching pursuit (MP) method, but better in that at each iteration, theresidual is recomputed using an orthogonal projection on the space of thepreviously chosen dictionary elements.

Examples

• Orthogonal Matching Pursuit

1.1.10.Bayesian Regression#

Bayesian regression techniques can be used to include regularizationparameters in the estimation procedure: the regularization parameter isnot set in a hard sense but tuned to the data at hand.

This can be done by introducinguninformative priorsover the hyper parameters of the model.The\(\ell_{2}\) regularization used inRidge regression and classification isequivalent to finding a maximum a posteriori estimation under a Gaussian priorover the coefficients\(w\) with precision\(\lambda^{-1}\).Instead of settinglambda manually, it is possible to treat it as a randomvariable to be estimated from the data.

To obtain a fully probabilistic model, the output\(y\) is assumedto be Gaussian distributed around\(X w\):

where\(\alpha\) is again treated as a random variable that is to beestimated from the data.

The advantages of Bayesian Regression are:

• It adapts to the data at hand.

• It can be used to include regularization parameters in theestimation procedure.

The disadvantages of Bayesian regression include:

• Inference of the model can be time consuming.

>>>fromsklearnimportlinear_model>>>X=[[0.,0.],[1.,1.],[2.,2.],[3.,3.]]>>>Y=[0.,1.,2.,3.]>>>reg=linear_model.BayesianRidge()>>>reg.fit(X,Y)BayesianRidge()

>>>reg.predict([[1,0.]])array([0.50000013])

>>>reg.coef_array([0.49999993, 0.49999993])

1.1.11.Logistic regression#

The logistic regression is implemented inLogisticRegression. Despiteits name, it is implemented as a linear model for classification rather thanregression in terms of the scikit-learn/ML nomenclature. The logisticregression is also known in the literature as logit regression,maximum-entropy classification (MaxEnt) or the log-linear classifier. In thismodel, the probabilities describing the possible outcomes of a single trialare modeled using alogistic function.

This implementation can fit binary, One-vs-Rest, or multinomial logisticregression with optional\(\ell_1\),\(\ell_2\) or Elastic-Netregularization.

Examples

• L1 Penalty and Sparsity in Logistic Regression

• Regularization path of L1- Logistic Regression

• Decision Boundaries of Multinomial and One-vs-Rest Logistic Regression

• Multiclass sparse logistic regression on 20newgroups

• MNIST classification using multinomial logistic + L1

• Plot classification probability

1.1.12.Generalized Linear Models#

Generalized Linear Models (GLM) extend linear models in two ways[10]. First, the predicted values\(\hat{y}\) are linked to a linearcombination of the input variables\(X\) via an inverse link function\(h\) as

Secondly, the squared loss function is replaced by the unit deviance\(d\) of a distribution in the exponential family (or more precisely, areproductive exponential dispersion model (EDM)[11]).

The minimization problem becomes:

where\(\alpha\) is the L2 regularization penalty. When sample weights areprovided, the average becomes a weighted average.

The following table lists some specific EDMs and their unit deviance :

The Probability Density Functions (PDF) of these distributions are illustratedin the following figure,

The Bernoulli distribution is a discrete probability distribution modelling aBernoulli trial - an event that has only two mutually exclusive outcomes.The Categorical distribution is a generalization of the Bernoulli distributionfor a categorical random variable. While a random variable in a Bernoullidistribution has two possible outcomes, a Categorical random variable can takeon one of K possible categories, with the probability of each categoryspecified separately.

The choice of the distribution depends on the problem at hand:

• If the target values\(y\) are counts (non-negative integer valued) orrelative frequencies (non-negative), you might use a Poisson distributionwith a log-link.

• If the target values are positive valued and skewed, you might try a Gammadistribution with a log-link.

• If the target values seem to be heavier tailed than a Gamma distribution, youmight try an Inverse Gaussian distribution (or even higher variance powers ofthe Tweedie family).

• If the target values\(y\) are probabilities, you can use the Bernoullidistribution. The Bernoulli distribution with a logit link can be used forbinary classification. The Categorical distribution with a softmax link can beused for multiclass classification.

References

McCullagh, Peter; Nelder, John (1989). Generalized Linear Models,Second Edition. Boca Raton: Chapman and Hall/CRC. ISBN 0-412-31760-5.

Jørgensen, B. (1992). The theory of exponential dispersion modelsand analysis of deviance. Monografias de matemática, no. 51. See alsoExponential dispersion model.

>>>fromsklearn.linear_modelimportTweedieRegressor>>>reg=TweedieRegressor(power=1,alpha=0.5,link='log')>>>reg.fit([[0,0],[0,1],[2,2]],[0,1,2])TweedieRegressor(alpha=0.5, link='log', power=1)>>>reg.coef_array([0.2463, 0.4337])>>>reg.intercept_np.float64(-0.7638)

1.1.13.Stochastic Gradient Descent - SGD#

Stochastic gradient descent is a simple yet very efficient approachto fit linear models. It is particularly useful when the number of samples(and the number of features) is very large.Thepartial_fit method allows online/out-of-core learning.

The classesSGDClassifier andSGDRegressor providefunctionality to fit linear models for classification and regressionusing different (convex) loss functions and different penalties.E.g., withloss="log",SGDClassifierfits a logistic regression model,while withloss="hinge" it fits a linear support vector machine (SVM).

You can refer to the dedicatedStochastic Gradient Descent documentation section for more details.

1.1.14.Perceptron#

ThePerceptron is another simple classification algorithm suitable forlarge scale learning. By default:

• It does not require a learning rate.

• It is not regularized (penalized).

• It updates its model only on mistakes.

The last characteristic implies that the Perceptron is slightly faster totrain than SGD with the hinge loss and that the resulting models aresparser.

In fact, thePerceptron is a wrapper around theSGDClassifierclass using a perceptron loss and a constant learning rate. Refer tomathematical section of the SGD procedurefor more details.

1.1.15.Passive Aggressive Algorithms#

The passive-aggressive algorithms are a family of algorithms for large-scalelearning. They are similar to the Perceptron in that they do not require alearning rate. However, contrary to the Perceptron, they include aregularization parameterC.

For classification,PassiveAggressiveClassifier can be used withloss='hinge' (PA-I) orloss='squared_hinge' (PA-II). For regression,PassiveAggressiveRegressor can be used withloss='epsilon_insensitive' (PA-I) orloss='squared_epsilon_insensitive' (PA-II).

1.1.16.Robustness regression: outliers and modeling errors#

Robust regression aims to fit a regression model in thepresence of corrupt data: either outliers, or error in the model.

1.1.16.2.RANSAC: RANdom SAmple Consensus#

RANSAC (RANdom SAmple Consensus) fits a model from random subsets ofinliers from the complete data set.

RANSAC is a non-deterministic algorithm producing only a reasonable result witha certain probability, which is dependent on the number of iterations (seemax_trials parameter). It is typically used for linear and non-linearregression problems and is especially popular in the field of photogrammetriccomputer vision.

The algorithm splits the complete input sample data into a set of inliers,which may be subject to noise, and outliers, which are e.g. caused by erroneousmeasurements or invalid hypotheses about the data. The resulting model is thenestimated only from the determined inliers.

Examples

• Robust linear model estimation using RANSAC

• Robust linear estimator fitting

Details of the algorithm#

Each iteration performs the following steps:

1. Selectmin_samples random samples from the original data and checkwhether the set of data is valid (seeis_data_valid).

2. Fit a model to the random subset (estimator.fit) and checkwhether the estimated model is valid (seeis_model_valid).

3. Classify all data as inliers or outliers by calculating the residualsto the estimated model (estimator.predict(X)-y) - all datasamples with absolute residuals smaller than or equal to theresidual_threshold are considered as inliers.

4. Save fitted model as best model if number of inlier samples ismaximal. In case the current estimated model has the same number ofinliers, it is only considered as the best model if it has better score.

These steps are performed either a maximum number of times (max_trials) oruntil one of the special stop criteria are met (seestop_n_inliers andstop_score). The final model is estimated using all inlier samples (consensusset) of the previously determined best model.

Theis_data_valid andis_model_valid functions allow to identify and rejectdegenerate combinations of random sub-samples. If the estimated model is notneeded for identifying degenerate cases,is_data_valid should be used as itis called prior to fitting the model and thus leading to better computationalperformance.

References#

• https://en.wikipedia.org/wiki/RANSAC

• “Random Sample Consensus: A Paradigm for Model Fitting with Applications toImage Analysis and Automated Cartography”Martin A. Fischler and Robert C. Bolles - SRI International (1981)

• “Performance Evaluation of RANSAC Family”Sunglok Choi, Taemin Kim and Wonpil Yu - BMVC (2009)

1.1.16.3.Theil-Sen estimator: generalized-median-based estimator#

TheTheilSenRegressor estimator uses a generalization of the median inmultiple dimensions. It is thus robust to multivariate outliers. Note howeverthat the robustness of the estimator decreases quickly with the dimensionalityof the problem. It loses its robustness properties and becomes nobetter than an ordinary least squares in high dimension.

Examples

• Theil-Sen Regression

• Robust linear estimator fitting

Theoretical considerations#

TheilSenRegressor is comparable to theOrdinary Least Squares(OLS) in terms of asymptotic efficiency and as anunbiased estimator. In contrast to OLS, Theil-Sen is a non-parametricmethod which means it makes no assumption about the underlyingdistribution of the data. Since Theil-Sen is a median-based estimator, itis more robust against corrupted data aka outliers. In univariatesetting, Theil-Sen has a breakdown point of about 29.3% in case of asimple linear regression which means that it can tolerate arbitrarycorrupted data of up to 29.3%.

The implementation ofTheilSenRegressor in scikit-learn follows ageneralization to a multivariate linear regression model[14] using thespatial median which is a generalization of the median to multipledimensions[15].

In terms of time and space complexity, Theil-Sen scales according to

\[\binom{n_{\text{samples}}}{n_{\text{subsamples}}}\]

which makes it infeasible to be applied exhaustively to problems with alarge number of samples and features. Therefore, the magnitude of asubpopulation can be chosen to limit the time and space complexity byconsidering only a random subset of all possible combinations.

References

[14]

Xin Dang, Hanxiang Peng, Xueqin Wang and Heping Zhang:Theil-Sen Estimators in a Multiple Linear Regression Model.

[15]

20. Kärkkäinen and S. Äyrämö:On Computation of Spatial Median for Robust Data Mining.

Also see theWikipedia page

1.1.16.4.Huber Regression#

TheHuberRegressor is different fromRidge because it applies alinear loss to samples that are defined as outliers by theepsilon parameter.A sample is classified as an inlier if the absolute error of that sample isless than the thresholdepsilon. It differs fromTheilSenRegressorandRANSACRegressor because it does not ignore the effect of the outliersbut gives a lesser weight to them.

Examples

• HuberRegressor vs Ridge on dataset with strong outliers

Mathematical details#

HuberRegressor minimizes

\[\min_{w, \sigma} {\sum_{i=1}^n\left(\sigma + H_{\epsilon}\left(\frac{X_{i}w - y_{i}}{\sigma}\right)\sigma\right) + \alpha {||w||_2}^2}\]

where the loss function is given by

\[\begin{split}H_{\epsilon}(z) = \begin{cases} z^2, & \text {if } |z| < \epsilon, \\ 2\epsilon|z| - \epsilon^2, & \text{otherwise}\end{cases}\end{split}\]

It is advised to set the parameterepsilon to 1.35 to achieve 95%statistical efficiency.

References

• Peter J. Huber, Elvezio M. Ronchetti: Robust Statistics, Concomitant scaleestimates, p. 172.

TheHuberRegressor differs from usingSGDRegressor with loss set tohuberin the following ways.

• HuberRegressor is scaling invariant. Onceepsilon is set, scalingX andydown or up by different values would produce the same robustness to outliers as before.as compared toSGDRegressor whereepsilon has to be set again whenX andy arescaled.

• HuberRegressor should be more efficient to use on data with small number ofsamples whileSGDRegressor needs a number of passes on the training data toproduce the same robustness.

Note that this estimator is different from theR implementation of RobustRegression because the Rimplementation does a weighted least squares implementation with weights given to eachsample on the basis of how much the residual is greater than a certain threshold.

1.1.17.Quantile Regression#

Quantile regression estimates the median or other quantiles of\(y\)conditional on\(X\), while ordinary least squares (OLS) estimates theconditional mean.

Quantile regression may be useful if one is interested in predicting aninterval instead of point prediction. Sometimes, prediction intervals arecalculated based on the assumption that prediction error is distributednormally with zero mean and constant variance. Quantile regression providessensible prediction intervals even for errors with non-constant (butpredictable) variance or non-normal distribution.

Based on minimizing the pinball loss, conditional quantiles can also beestimated by models other than linear models. For example,GradientBoostingRegressor can predict conditionalquantiles if its parameterloss is set to"quantile" and parameteralpha is set to the quantile that should be predicted. See the example inPrediction Intervals for Gradient Boosting Regression.

Most implementations of quantile regression are based on linear programmingproblem. The current implementation is based onscipy.optimize.linprog.

Examples

• Quantile regression

1.1.18.Polynomial regression: extending linear models with basis functions#

One common pattern within machine learning is to use linear models trainedon nonlinear functions of the data. This approach maintains the generallyfast performance of linear methods, while allowing them to fit a much widerrange of data.

Here is an example of applying this idea to one-dimensional data, usingpolynomial features of varying degrees:

This figure is created using thePolynomialFeatures transformer, whichtransforms an input data matrix into a new data matrix of a given degree.It can be used as follows:

>>>fromsklearn.preprocessingimportPolynomialFeatures>>>importnumpyasnp>>>X=np.arange(6).reshape(3,2)>>>Xarray([[0, 1],       [2, 3],       [4, 5]])>>>poly=PolynomialFeatures(degree=2)>>>poly.fit_transform(X)array([[ 1.,  0.,  1.,  0.,  0.,  1.],       [ 1.,  2.,  3.,  4.,  6.,  9.],       [ 1.,  4.,  5., 16., 20., 25.]])

The features ofX have been transformed from\([x_1, x_2]\) to\([1, x_1, x_2, x_1^2, x_1 x_2, x_2^2]\), and can now be used withinany linear model.

This sort of preprocessing can be streamlined with thePipeline tools. A single object representing a simplepolynomial regression can be created and used as follows:

>>>fromsklearn.preprocessingimportPolynomialFeatures>>>fromsklearn.linear_modelimportLinearRegression>>>fromsklearn.pipelineimportPipeline>>>importnumpyasnp>>>model=Pipeline([('poly',PolynomialFeatures(degree=3)),...('linear',LinearRegression(fit_intercept=False))])>>># fit to an order-3 polynomial data>>>x=np.arange(5)>>>y=3-2*x+x**2-x**3>>>model=model.fit(x[:,np.newaxis],y)>>>model.named_steps['linear'].coef_array([ 3., -2.,  1., -1.])

The linear model trained on polynomial features is able to exactly recoverthe input polynomial coefficients.

In some cases it’s not necessary to include higher powers of any single feature,but only the so-calledinteraction featuresthat multiply together at most\(d\) distinct features.These can be gotten fromPolynomialFeatures with the settinginteraction_only=True.

For example, when dealing with boolean features,\(x_i^n = x_i\) for all\(n\) and is therefore useless;but\(x_i x_j\) represents the conjunction of two booleans.This way, we can solve the XOR problem with a linear classifier:

>>>fromsklearn.linear_modelimportPerceptron>>>fromsklearn.preprocessingimportPolynomialFeatures>>>importnumpyasnp>>>X=np.array([[0,0],[0,1],[1,0],[1,1]])>>>y=X[:,0]^X[:,1]>>>yarray([0, 1, 1, 0])>>>X=PolynomialFeatures(interaction_only=True).fit_transform(X).astype(int)>>>Xarray([[1, 0, 0, 0],       [1, 0, 1, 0],       [1, 1, 0, 0],       [1, 1, 1, 1]])>>>clf=Perceptron(fit_intercept=False,max_iter=10,tol=None,...shuffle=False).fit(X,y)

And the classifier “predictions” are perfect:

>>>clf.predict(X)array([0, 1, 1, 0])>>>clf.score(X,y)1.0