Thank you for all the pointers @lorentzenchr !

We need to update test_regression to run a random forest with the poisson criterion. The target would need to be adjust so that it is greater than zero.

Do you mean similar to:

• scikit-learn/sklearn/tree/tests/test_tree.py

  Line 2026 in ab65c8b

  def test_poisson_vs_mse():

• scikit-learn/sklearn/tree/tests/test_tree.py

  Line 1980 in ab65c8b

  def test_balance_property(criterion, Tree):

What about

• scikit-learn/sklearn/tree/tests/test_tree.py

  Line 1993 in ab65c8b

  def test_poisson_zero_nodes(seed):

I am thinking of doing simple integration tests, such as:

assert np.all(reg.predict(X) > 0)

At a glance, due to the bagging nature of forests, I do not think all the conditions that held for single trees also holds for forests.

test_poisson_vs_mse could work, but I suspect some fine-tuning.

apologies, I'm not quite following the convo - for now, I've tried to rework check_regression_criterion so that it can run with the Poisson criterion.

def check_regression_criterion(name, criterion):
    # Check consistency on regression dataset.
    ForestRegressor = FOREST_REGRESSORS[name]

    if criterion == 'poisson':
        mask = y_reg > 0
        y_regr = y_reg * mask + 1
    else:
        y_regr = y_reg

    reg = ForestRegressor(n_estimators=5, criterion=criterion,
                          random_state=1)
    reg.fit(X_reg, y_regr)
    score = reg.score(X_reg, y_regr)
    assert score > 0.93, ("Failed with max_features=None, criterion %s "
                          "and score = %f" % (criterion, score))

    reg = ForestRegressor(n_estimators=5, criterion=criterion,
                          max_features=6, random_state=1)
    reg.fit(X_reg, y_regr)
    score = reg.score(X_reg, y_regr)
    assert score > 0.92, ("Failed with max_features=6, criterion %s "
                          "and score = %f" % (criterion, score))

This might not be what we're looking for, but thought I'd put something out to try and improve my understanding.

@bsun94 The point is that we should test, on some dataset, that the poisson forest regression actually works. One possibility is to test that a poisson forest has better predictive performance on a poisson distributed target compared to as squared error forest, see test_poisson_vs_mse in the tree module. (Note that this test required some fine tuning for single trees.)
Another possibility is to check the balance property, i.e. sum(y_true) == sum(y_pred) on the training set. This should also be true for squared error.

Ideally, we have both tests.

I see, thanks Christian. I've replicated the test_poisson_vs_mse test in the test_forest module - adapted to use RandomForestRegressor. I've tried playing around with n_features, different settings in rng.uniform and rng.poisson, but it seems no matter what I try, the test fails at metric_poi < metric_mse - mse having the lower error. Are there other ways for me to fine-tune/am I going about this the right way?

Taking test_poisson_vs_mse and replacing the tree models by

forest_poi = RandomForestRegressor(
    criterion="poisson",
    n_estimators=10,
    bootstrap=False,
    min_samples_split=10,
    random_state=rng
)
forest_mse = RandomForestRegressor(
    criterion="mse",
    n_estimators=10,
    bootstrap=False,
    min_samples_split=10,
    random_state=rng
)

works for me. See details.

import numpy as np
from sklearn.dummy import DummyRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import mean_squared_error
from sklearn.metrics import mean_poisson_deviance
from sklearn.model_selection import train_test_split
from sklearn.tree import DecisionTreeRegressor
from sklearn import datasets


rng = np.random.RandomState(42)
n_train, n_test, n_features = 500, 500, 10
X = datasets.make_low_rank_matrix(n_samples=n_train + n_test,
                                  n_features=n_features, random_state=rng)
# We create a log-linear Poisson model and downscale coef as it will get
# exponentiated.
coef = rng.uniform(low=-2, high=2, size=n_features) / np.max(X, axis=0)
y = rng.poisson(lam=np.exp(X @ coef))
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=n_test,
                                                    random_state=rng)

forest_poi = RandomForestRegressor(
    criterion="poisson",
    n_estimators=10,
    bootstrap=False,
    random_state=rng
)
forest_mse = RandomForestRegressor(
    criterion="mse",
    n_estimators=10,
    bootstrap=False,
    random_state=rng
)
forest_poi.fit(X_train, y_train)
forest_mse.fit(X_train, y_train)
dummy = DummyRegressor(strategy="mean").fit(X_train, y_train)

for X, y, val in [(X_train, y_train, "train"), (X_test, y_test, "test")]:
    metric_poi = mean_poisson_deviance(y, forest_poi.predict(X))
    # squared_error might produce non-positive predictions => clip
    metric_mse = mean_poisson_deviance(y, np.clip(forest_mse.predict(X),
                                                  1e-15, None))
    metric_dummy = mean_poisson_deviance(y, dummy.predict(X))
    # As squared_error might correctly predict 0 in train set, its train
    # score can be better than Poisson. This is no longer the case for the
    # test set.
    if val == "test":
        assert metric_poi < metric_mse
    assert metric_poi < metric_dummy

No worries @lorentzenchr - thank you for your guidance throughout!
To be honest, I don't think I understand all the Maths here haha (why does this work?) - would there be any online resources you could refer me to for me to read up?

@bsun94 Which aspect to you mean? The loss functions aka splitting criteria, how single decision trees work or how the random forest makes use of many single trees (with nice theoretical properties)?

@bsun94 Which aspect to you mean? The loss functions aka splitting criteria, how single decision trees work or how the random forest makes use of many single trees (with nice theoretical properties)?

@lorentzenchr On all three would be great haha

@bsun94 To advertise this nice library, why don't you have a look at the user guide:

• Decision Trees: https://scikit-learn.org/stable/modules/tree.html
• Random Forests: https://scikit-learn.org/stable/modules/ensemble.html#forests-of-randomized-trees

I don't have a good reference for loss functions without knowing your background. In principle, it has components from probability theory (Poisson distribution), statistics (estimation theory, Poisson GLM), statistical learning theory and forecasting (and maybe more).

@bsun94 Could you resolve merge conflicts by merging the main branch?

thank you @lorentzenchr and @thomasjpfan for all your guidance in this!

@bsun94 Happy to help out. Maybe there is a next time:smirk:

I don't see why not!