Counting on Poisson Regression with

Thomas J. Fan
@thomasjpfan
This talk on Github: thomasjpfan/scipy-2022-poisson

Agriculture

Risk modeling

Predictive maintenance

Poisson Regression!

PoissonRegressor

from sklearn.linear_model import PoissonRegressor
reg = PoissonRegressor()

Poisson Regression!

PoissonRegressor

from sklearn.linear_model import PoissonRegressor
reg = PoissonRegressor()

RandomForestRegressor

from sklearn.ensemble import RandomForestRegressor
reg = RandomForestRegressor(criterion="poisson")

Poisson Regression!

PoissonRegressor

from sklearn.linear_model import PoissonRegressor
reg = PoissonRegressor()

RandomForestRegressor

from sklearn.ensemble import RandomForestRegressor
reg = RandomForestRegressor(criterion="poisson")

HistGradientBoostingRegressor

from sklearn.ensemble import HistGradientBoostingRegressor
reg = HistGradientBoostingRegressor(loss="poisson")

ColumnTransformer

preprocessor = ColumnTransformer([
    (
        "cyclic_hour",
        SplineTransformer(n_knots=13, extrapolation="periodic"),
        ["hour"]
    ),
    (
        "categorical",
        OneHotEncoder(handle_unknown="ignore"),
        ["is_holiday", "weather_code", "is_weekend", "season"]
    ),
], remainder=MinMaxScaler())

ColumnTransformer

preprocessor = ColumnTransformer([
    (
        "cyclic_hour",
        SplineTransformer(n_knots=13, extrapolation="periodic"),
        ["hour"]
    ),
    (
        "categorical",
        OneHotEncoder(handle_unknown="ignore"),
        ["is_holiday", "weather_code", "is_weekend", "season"]
    ),
], remainder=MinMaxScaler())

ColumnTransformer

preprocessor = ColumnTransformer([
    (
        "cyclic_hour",
        SplineTransformer(n_knots=13, extrapolation="periodic"),
        ["hour"]
    ),
    (
        "categorical",
        OneHotEncoder(handle_unknown="ignore"),
        ["is_holiday", "weather_code", "is_weekend", "season"]
    ),
], remainder=MinMaxScaler())

ColumnTransformer

preprocessor = ColumnTransformer([
    (
        "cyclic_hour",
        SplineTransformer(n_knots=13, extrapolation="periodic"),
        ["hour"]
    ),
    (
        "categorical",
        OneHotEncoder(handle_unknown="ignore"),
        ["is_holiday", "weather_code", "is_weekend", "season"]
    ),
], remainder=MinMaxScaler())

Periodic spline features

Time-related feature engineering example

PoissonRegressor 🎢

Generalized Linear Models (GLM)

$\hat{y}(w, X) = h(Xw)$

where $\hat{y}$ is the predicted values, $X$ are features, and $h$ is the inverse link function.

Generalized Linear Models (GLM)

$\hat{y}(w, X) = h(Xw)$

where $\hat{y}$ is the predicted values, $X$ are features, and $h$ is the inverse link function.

Minimization problem becomes:

$\min_{w}\frac{1}{2n} \sum_i d(y_i, \hat{y}_i) + \frac{\alpha}{2}||w||_2^2$

where $\alpha$ is the L2 regularization penalty.

User Guide

Deviance

User Guide

Minimization Problem

scipy.optimize.minimize is used with L-BFGS-B

$\min_{w}\frac{1}{2n} \sum_i d(y_i, \hat{y}_i) + \frac{\alpha}{2}||w||_2^2$

Minimization Problem

scipy.optimize.minimize is used with L-BFGS-B

$\min_{w}\frac{1}{2n} \sum_i d(y_i, \hat{y}_i) + \frac{\alpha}{2}||w||_2^2$

Cholesky based Newton solver: PR #23314
Newton-LSMR: PR #23507

Minimization Problem

scipy.optimize.minimize is used with L-BFGS-B

$\min_{w}\frac{1}{2n} \sum_i d(y_i, \hat{y}_i) + \frac{\alpha}{2}||w||_2^2$

Cholesky based Newton solver: PR #23314
Newton-LSMR: PR #23507

Regularization by Default!

PoissonRegressor(alpha=1.0)

Preprocessor + Linear Model

Poisson

poisson = make_pipeline(preprocessor, PoissonRegressor(...))

Preprocessor + Linear Model

Poisson

poisson = make_pipeline(preprocessor, PoissonRegressor(...))

Ridge

ridge = make_pipeline(preprocessor, Ridge())

Evaluation

from sklearn.model_selection import TimeSeriesSplit
cv = TimeSeriesSplit(
    n_splits=50,
    max_train_size=10000,
    test_size=336,
)

Results - Linear Models

Distributions - Linear Models

Calibration for Regression

Calibration - Linear Models

Random Forest 🎄🎄🎄

Random Forest With Poisson 🎄🎄🎄

RandomForestRegressor(criterion="poisson")

How does Poisson Influence the Random Forest?

$\text{H}(Q) = \dfrac{1}{n}\sum_{i\in Q} y_i * \log\left(\dfrac{y_i}{\hat{y}_i}\right) + \hat{y}_i - y_i$

where

$y_i$ is the true value
$\hat{y}_i$ is the predicted value
$Q$ is a node,
$n$ is the number of data points in node.

Details in User Guide

Random Forest

Results - Random Forest

Distributions - Random Forest

Calibration - Random Forest

Histogram-based Gradient Boosting Trees 🏂

HistGradientBoostingRegressor With Poisson 🏂

HistGradientBoostingRegressor(loss="poisson")

HistGradientBoostingRegressor Overview 🏂

Initial Condition

$\hat{f}^{(0)} = C$

where $C$ is the maximum likelihood estimate for loss.

HistGradientBoostingRegressor Overview 🏂

Initial Condition

$\hat{f}^{(0)} = C$

where $C$ is the maximum likelihood estimate for loss.

Iterations

$\hat{f}^{(t)}=\hat{f}^{(t-1)} + \nu \hat{h}^{(t)}$

where

$\hat{f}^{(t)}$ is the regressor at iteration $t$
$\nu$ are learning rate
$\hat{h}^{(t)}$ are trees using the gradient and hessians.

How does Poisson Influence the Algorithm?

Growing Trees $\hat{h}^{(t)}$ by to evaluating splits:

$\text{Gain} = \frac{1}{2}\left[\dfrac{G_L^2}{H_L+\lambda} + \dfrac{G_R^2}{H_R+\lambda} - \frac{(G_L+G_R)^2}{H_L+H_R+\lambda}\right]$

where

$\lambda$ is the l2 regularization
$G_L$ and $G_R$ is sum of the gradients
$H_L$ and $H_R$ is sum of the hessians

How does Poisson Influence the Algorithm?

Growing Trees $\hat{h}^{(t)}$ by to evaluating splits:

$\text{Gain} = \frac{1}{2}\left[\dfrac{G_L^2}{H_L+\lambda} + \dfrac{G_R^2}{H_R+\lambda} - \frac{(G_L+G_R)^2}{H_L+H_R+\lambda}\right]$

where

$\lambda$ is the l2 regularization
$G_L$ and $G_R$ is sum of the gradients
$H_L$ and $H_R$ is sum of the hessians

More Details @

SciPy 2019 Fast Gradient Boosting Decision Trees with PyGBM and Numba

Linking $f$ with $y$

$\hat{y}^{(t)} = h(\hat{f}^{(t)})$

where $h$ is the inverse link function.

Linking $f$ with $y$

$\hat{y}^{(t)} = h(\hat{f}^{(t)})$

where $h$ is the inverse link function.

Poisson's Inverse Link function

$h(z) = \exp(z)$

Linking $f$ with $y$

$\hat{y}^{(t)} = h(\hat{f}^{(t)})$

where $h$ is the inverse link function.

Poisson's Inverse Link function

$h(z) = \exp(z)$

Looks link the GLMs

$\hat{y}(w, X) = h(Xw)$

Results - Hist Gradient Boosting

Results - Hist Gradient Boosting 🔎

Distributions - Hist Gradient Boosting

Calibration - Hist Gradient Boosting

Example of Predictions

Poisson Regression with Bike Share Data 🚲🚲🚲PoissonRegressor()RandomForestRegressor(criterion="poisson")HistGradientBoostingRegressor(loss="poisson")

Two More Topics 🔎

Zero-Inflated Poisson Regression

Zero-Inflated Poisson RegressionScikit-lego!from sklego.meta import ZeroInflatedRegressor
poisson_zero =  ZeroInflatedRegressor(
     classifier=HistGradientBoostingClassifier(),
     regressor=HistGradientBoostingRegressor(loss="poisson"),
)

Exposure

df["Frequency"] = df["ClaimNb"] / df["Exposure"]
poisson_gbrt.fit(
    df_train, df_train["Frequency"],
    regressor__sample_weight=df_train["Exposure"],
)

Poisson regression and non-normal loss example

Counting on Poisson Regression

`PoissonRegressor()`

`RandomForestRegressor(criterion="poisson")`

`HistGradientBoostingRegressor(loss="poisson")`

Thomas J. Fan
@thomasjpfan

This workshop on Github: github.com/thomasjpfan/scipy-2022-poisson

Help

Keyboard shortcuts

↑, ←, Pg Up, k

Go to previous slide

↓, →, Pg Dn, Space, j

Go to next slide

Home

Go to first slide

End

Go to last slide

Number + Return

Go to specific slide

b / m / f

Toggle blackout / mirrored / fullscreen mode

Clone slideshow

Toggle presenter mode

Restart the presentation timer

?, h

Toggle this help

Counting on Poisson Regression with

Agriculture

Risk modeling

Predictive maintenance

Poisson Regression!

PoissonRegressor

Poisson Regression!

PoissonRegressor

RandomForestRegressor

Poisson Regression!

PoissonRegressor

RandomForestRegressor

HistGradientBoostingRegressor

Bike Sharing Dataset 🚲

Bike Sharing Dataset 🚲

Bike Sharing Dataset 🚲

Bike Sharing Dataset 🚲

ColumnTransformer

ColumnTransformer

ColumnTransformer

ColumnTransformer

Periodic spline features

PoissonRegressor 🎢

Generalized Linear Models (GLM)

Generalized Linear Models (GLM)

Minimization problem becomes:

Deviance

Minimization Problem

Minimization Problem

Minimization Problem

Regularization by Default!

Preprocessor + Linear Model

Poisson

Preprocessor + Linear Model

Poisson

Ridge

Evaluation

Results - Linear Models

Distributions - Linear Models

Calibration for Regression

Calibration for Regression

Calibration for Regression

Calibration for Regression

Calibration - Linear Models

Random Forest 🎄🎄🎄

Random Forest With Poisson 🎄🎄🎄

How does Poisson Influence the Random Forest?

Random Forest

Results - Random Forest

Distributions - Random Forest

Calibration - Random Forest

Histogram-based Gradient Boosting Trees 🏂

HistGradientBoostingRegressor With Poisson 🏂

HistGradientBoostingRegressor Overview 🏂

Initial Condition

HistGradientBoostingRegressor Overview 🏂

Initial Condition

Iterations

How does Poisson Influence the Algorithm?

Growing Trees h^(t)\hat{h}^{(t)}h^(t) by to evaluating splits:

How does Poisson Influence the Algorithm?

Growing Trees h^(t)\hat{h}^{(t)}h^(t) by to evaluating splits:

More Details @

Linking fff with yyy

Linking fff with yyy

Poisson's Inverse Link function

Linking fff with yyy

Poisson's Inverse Link function

Looks link the GLMs

Results - Hist Gradient Boosting

Results - Hist Gradient Boosting 🔎

Distributions - Hist Gradient Boosting

Calibration - Hist Gradient Boosting

Example of Predictions

Poisson Regression with Bike Share Data 🚲🚲🚲

PoissonRegressor()

RandomForestRegressor(criterion="poisson")

HistGradientBoostingRegressor(loss="poisson")

Two More Topics 🔎

Zero-Inflated Poisson Regression

Growing Trees $\hat{h}^{(t)}$ by to evaluating splits:

Growing Trees $\hat{h}^{(t)}$ by to evaluating splits:

Linking $f$ with $y$

Linking $f$ with $y$

Linking $f$ with $y$

`PoissonRegressor()`

`RandomForestRegressor(criterion="poisson")`

`HistGradientBoostingRegressor(loss="poisson")`

`PoissonRegressor()`

`RandomForestRegressor(criterion="poisson")`

`HistGradientBoostingRegressor(loss="poisson")`