added "quantile" family

Author: mathias-von-ottenbreit

API_REFERENCE.md

@@ -1,6 +1,6 @@
 # APLRRegressor
 
-## class aplr.APLRRegressor(m:int=1000, v:float=0.1, random_state:int=0, family:str="gaussian", link_function:str="identity", n_jobs:int=0, validation_ratio:float=0.2, intercept:float=np.nan, bins:int=300, max_interaction_level:int=1, max_interactions:int=100000, min_observations_in_split:int=20, ineligible_boosting_steps_added:int=10, max_eligible_terms:int=5, verbosity:int=0, tweedie_power:float=1.5, validation_tuning_metric:str="default")
+## class aplr.APLRRegressor(m:int=1000, v:float=0.1, random_state:int=0, family:str="gaussian", link_function:str="identity", n_jobs:int=0, validation_ratio:float=0.2, intercept:float=np.nan, bins:int=300, max_interaction_level:int=1, max_interactions:int=100000, min_observations_in_split:int=20, ineligible_boosting_steps_added:int=10, max_eligible_terms:int=5, verbosity:int=0, tweedie_power:float=1.5, validation_tuning_metric:str="default", quantile:float=0.5)
 
 ### Constructor parameters
 
@@ -14,7 +14,7 @@ The learning rate. Must be greater than zero and not more than one. The higher t
 Used to randomly split training observations into training and validation if ***validation_set_indexes*** is not specified when fitting.
 
 #### family (default = "gaussian")
-Determines the loss function used. Allowed values are "gaussian", "binomial", "poisson", "gamma", "tweedie", "group_gaussian" and "mae". This is used together with ***link_function***. When ***family*** is "group_gaussian" then the "group" argument in the ***fit*** method must be provided. In the latter case APLR will try to minimize group MSE when training the model.
+Determines the loss function used. Allowed values are "gaussian", "binomial", "poisson", "gamma", "tweedie", "group_gaussian", "mae" and "quantile". This is used together with ***link_function***. When ***family*** is "group_gaussian" then the "group" argument in the ***fit*** method must be provided. In the latter case APLR will try to minimize group MSE when training the model. The ***family*** "quantile" is used together with the ***quantile*** constructor parameter.
 
 #### link_function (default = "identity")
 Determines how the linear predictor is transformed to predictions. Allowed values are "identity", "logit" and "log". For an ordinary regression model use ***family*** "gaussian" and ***link_function*** "identity". For logistic regression use ***family*** "binomial" and ***link_function*** "logit". For a multiplicative model use the "log" ***link_function***. The "log" ***link_function*** often works best with a "poisson", "gamma" or "tweedie" ***family***, depending on the data. The ***family*** "poisson", "gamma" or "tweedie" should only be used with the "log" ***link_function***. Inappropriate combinations of ***family*** and ***link_function*** may result in a warning message when fitting the model and/or a poor model fit. Please note that values other than "identity" typically require a significantly higher ***m*** (or ***v***) in order to converge.
@@ -55,6 +55,10 @@ Specifies the variance power for the "tweedie" ***family***.
 #### validation_tuning_metric (default = "default")
 Specifies which metric to use for validating the model and tuning ***m***. Available options are "default" (using the same methodology as when calculating the training error), "mse", "mae", "negative_gini" and "rankability". The default is often a choice that fits well with respect to the ***family*** chosen. However, if you want to use ***family*** or ***tweedie_power*** as tuning parameters then the default is not suitable. "rankability" uses a methodology similar to the one described in https://towardsdatascience.com/how-to-calculate-roc-auc-score-for-regression-models-c0be4fdf76bb except that the metric is inverted and can be weighted by sample weights.
 
+#### quantile (default = 0.5)
+Specifies the quantile to use when ***family*** is "quantile".
+
+
 ## Method: fit(X:npt.ArrayLike, y:npt.ArrayLike, sample_weight:npt.ArrayLike = np.empty(0), X_names:List[str]=[], validation_set_indexes:List[int]=[], prioritized_predictors_indexes:List[int]=[], monotonic_constraints:List[int]=[], group:npt.ArrayLike = np.empty(0))
 
 ***This method fits the model to data.***

aplr/aplr.py

@@ -5,7 +5,7 @@
 
 
 class APLRRegressor():
-    def __init__(self, m:int=1000, v:float=0.1, random_state:int=0, family:str="gaussian", link_function:str="identity", n_jobs:int=0, validation_ratio:float=0.2, intercept:float=np.nan, bins:int=300, max_interaction_level:int=1, max_interactions:int=100000, min_observations_in_split:int=20, ineligible_boosting_steps_added:int=10, max_eligible_terms:int=5, verbosity:int=0, tweedie_power:float=1.5, validation_tuning_metric:str="default"):
+    def __init__(self, m:int=1000, v:float=0.1, random_state:int=0, family:str="gaussian", link_function:str="identity", n_jobs:int=0, validation_ratio:float=0.2, intercept:float=np.nan, bins:int=300, max_interaction_level:int=1, max_interactions:int=100000, min_observations_in_split:int=20, ineligible_boosting_steps_added:int=10, max_eligible_terms:int=5, verbosity:int=0, tweedie_power:float=1.5, validation_tuning_metric:str="default", quantile:float=0.5):
         self.m=m
         self.v=v
         self.random_state=random_state
@@ -23,6 +23,7 @@ def __init__(self, m:int=1000, v:float=0.1, random_state:int=0, family:str="gaus
         self.verbosity=verbosity
         self.tweedie_power=tweedie_power
         self.validation_tuning_metric=validation_tuning_metric
+        self.quantile=quantile
 
         #Creating aplr_cpp and setting parameters
         self.APLRRegressor=aplr_cpp.APLRRegressor()
@@ -47,6 +48,7 @@ def __set_params_cpp(self):
         self.APLRRegressor.verbosity=self.verbosity
         self.APLRRegressor.tweedie_power=self.tweedie_power
         self.APLRRegressor.validation_tuning_metric=self.validation_tuning_metric
+        self.APLRRegressor.quantile=self.quantile
 
     def fit(self, X:npt.ArrayLike, y:npt.ArrayLike, sample_weight:npt.ArrayLike = np.empty(0), X_names:List[str]=[], validation_set_indexes:List[int]=[], prioritized_predictors_indexes:List[int]=[], monotonic_constraints:List[int]=[],group:npt.ArrayLike = np.empty(0)):
         self.__set_params_cpp()
@@ -116,7 +118,8 @@ def get_params(self, deep=True):
             "ineligible_boosting_steps_added":self.ineligible_boosting_steps_added,
             "max_eligible_terms":self.max_eligible_terms,
             "tweedie_power":self.tweedie_power,
-            "validation_tuning_metric":self.validation_tuning_metric
+            "validation_tuning_metric":self.validation_tuning_metric,
+            "quantile":self.quantile
         }
 
     #For sklearn

cpp/APLRRegressor.h

@@ -141,13 +141,14 @@ class APLRRegressor
     double max_training_prediction_or_response;
     std::vector<size_t> validation_indexes;
     std::string validation_tuning_metric;
+    double quantile;
 
     //Methods
     APLRRegressor(size_t m=1000,double v=0.1,uint_fast32_t random_state=std::numeric_limits<uint_fast32_t>::lowest(),std::string family="gaussian",
         std::string link_function="identity", size_t n_jobs=0, double validation_ratio=0.2,double intercept=NAN_DOUBLE,
         size_t reserved_terms_times_num_x=100, size_t bins=300,size_t verbosity=0,size_t max_interaction_level=1,size_t max_interactions=100000,
         size_t min_observations_in_split=20, size_t ineligible_boosting_steps_added=10, size_t max_eligible_terms=5,double tweedie_power=1.5,
-        std::string validation_tuning_metric="default");
+        std::string validation_tuning_metric="default", double quantile=0.5);
     APLRRegressor(const APLRRegressor &other);
     ~APLRRegressor();
     void fit(const MatrixXd &X,const VectorXd &y,const VectorXd &sample_weight=VectorXd(0),const std::vector<std::string> &X_names={},const std::vector<size_t> &validation_set_indexes={},
@@ -174,15 +175,15 @@ class APLRRegressor
 APLRRegressor::APLRRegressor(size_t m,double v,uint_fast32_t random_state,std::string family,std::string link_function,size_t n_jobs,
     double validation_ratio,double intercept,size_t reserved_terms_times_num_x,size_t bins,size_t verbosity,size_t max_interaction_level,
     size_t max_interactions,size_t min_observations_in_split,size_t ineligible_boosting_steps_added,size_t max_eligible_terms,double tweedie_power,
-    std::string validation_tuning_metric):
+    std::string validation_tuning_metric, double quantile):
         reserved_terms_times_num_x{reserved_terms_times_num_x},intercept{intercept},m{m},v{v},
         family{family},link_function{link_function},validation_ratio{validation_ratio},n_jobs{n_jobs},random_state{random_state},
         bins{bins},verbosity{verbosity},max_interaction_level{max_interaction_level},
         intercept_steps{VectorXd(0)},max_interactions{max_interactions},interactions_eligible{0},validation_error_steps{VectorXd(0)},
         min_observations_in_split{min_observations_in_split},ineligible_boosting_steps_added{ineligible_boosting_steps_added},
         max_eligible_terms{max_eligible_terms},number_of_base_terms{0},tweedie_power{tweedie_power},min_training_prediction_or_response{NAN_DOUBLE},
         max_training_prediction_or_response{NAN_DOUBLE}, validation_tuning_metric{validation_tuning_metric},
-        validation_indexes{std::vector<size_t>(0)}
+        validation_indexes{std::vector<size_t>(0)}, quantile{quantile}
 {
 }
 
@@ -198,7 +199,7 @@ APLRRegressor::APLRRegressor(const APLRRegressor &other):
     max_eligible_terms{other.max_eligible_terms},number_of_base_terms{other.number_of_base_terms},
     feature_importance{other.feature_importance},tweedie_power{other.tweedie_power},min_training_prediction_or_response{other.min_training_prediction_or_response},
     max_training_prediction_or_response{other.max_training_prediction_or_response},validation_tuning_metric{other.validation_tuning_metric},
-    validation_indexes{other.validation_indexes}
+    validation_indexes{other.validation_indexes}, quantile{other.quantile}
 {
 }
 
@@ -250,6 +251,8 @@ void APLRRegressor::throw_error_if_family_does_not_exist()
         family_exists=true;
     else if(family=="mae")
         family_exists=true;
+    else if(family=="quantile")
+        family_exists=true;
     if(!family_exists)
         throw std::runtime_error("Family "+family+" is not available in APLR.");   
 }
@@ -613,6 +616,18 @@ VectorXd APLRRegressor::calculate_neg_gradient_current(const VectorXd &sample_we
         double mae{calculate_errors(y_train,predictions_current,sample_weight_train,"mae").mean()};
         output=(y_train.array() - predictions_current.array()).sign()*mae;
     }
+    else if(family=="quantile")
+    {
+        double mae{calculate_errors(y_train,predictions_current,sample_weight_train,"mae").mean()};
+        output=(y_train.array() - predictions_current.array()).sign()*mae;
+        for (Eigen::Index i = 0; i < y_train.size(); ++i)
+        {
+            if(y_train[i]<predictions_current[i])
+                output[i] *= 1-quantile;
+            else
+                output[i] *= quantile;
+        }
+    }    
 
     if(link_function!="identity")
         output=output.array()*differentiate_predictions().array();
@@ -1056,7 +1071,7 @@ void APLRRegressor::calculate_and_validate_validation_error(size_t boosting_step
 void APLRRegressor::calculate_validation_error(size_t boosting_step, const VectorXd &predictions)
 {
     if(validation_tuning_metric=="default")
-        validation_error_steps[boosting_step]=calculate_mean_error(calculate_errors(y_validation,predictions,sample_weight_validation,family,tweedie_power,group_validation,unique_groups_validation),sample_weight_validation);
+        validation_error_steps[boosting_step]=calculate_mean_error(calculate_errors(y_validation,predictions,sample_weight_validation,family,tweedie_power,group_validation,unique_groups_validation,quantile),sample_weight_validation);
     else if(validation_tuning_metric=="mse")
         validation_error_steps[boosting_step]=calculate_mean_error(calculate_errors(y_validation,predictions,sample_weight_validation,FAMILY_GAUSSIAN),sample_weight_validation);
     else if(validation_tuning_metric=="mae")

cpp/functions.h

@@ -133,7 +133,22 @@ VectorXd calculate_absolute_errors(const VectorXd &y,const VectorXd &predicted)
     return errors;
 }
 
-VectorXd calculate_errors(const VectorXd &y,const VectorXd &predicted,const VectorXd &sample_weight=VectorXd(0),const std::string &family="gaussian",double tweedie_power=1.5, const VectorXi &group=VectorXi(0), const std::set<int> &unique_groups={})
+VectorXd calculate_quantile_errors(const VectorXd &y,const VectorXd &predicted, double quantile)
+{
+    VectorXd errors{calculate_absolute_errors(y,predicted)};
+    for (Eigen::Index i = 0; i < y.size(); ++i)
+    {
+        if(y[i]<predicted[i])
+            errors[i] *= 1-quantile;
+        else
+            errors[i] *= quantile;
+    }
+    
+    return errors;
+}
+
+VectorXd calculate_errors(const VectorXd &y,const VectorXd &predicted,const VectorXd &sample_weight=VectorXd(0),const std::string &family="gaussian",
+    double tweedie_power=1.5, const VectorXi &group=VectorXi(0), const std::set<int> &unique_groups={}, double quantile=0.5)
 {   
     VectorXd errors;
     if(family=="gaussian")
@@ -150,6 +165,8 @@ VectorXd calculate_errors(const VectorXd &y,const VectorXd &predicted,const Vect
         errors=calculate_group_gaussian_errors(y,predicted,group,unique_groups);
     else if(family=="mae")
         errors=calculate_absolute_errors(y,predicted);
+    else if(family=="quantile")
+        errors=calculate_quantile_errors(y,predicted,quantile);
 
     if(sample_weight.size()>0)
         errors=errors.array()*sample_weight.array();

cpp/pythonbinding.cpp

@@ -11,14 +11,16 @@ namespace py = pybind11;
 
 PYBIND11_MODULE(aplr_cpp, m) {
     py::class_<APLRRegressor>(m, "APLRRegressor",py::module_local())
-        .def(py::init<int&, double&, int&, std::string&,std::string&,int&,double&,double&,int&,int&,int&,int&,int&,int&,int&,int&,double&,std::string&>(),
+        .def(py::init<int&, double&, int&, std::string&,std::string&,int&,double&,double&,int&,int&,int&,int&,int&,int&,int&,int&,double&,std::string&,
+            double&>(),
             py::arg("m")=1000,py::arg("v")=0.1,py::arg("random_state")=0,py::arg("family")="gaussian",py::arg("link_function")="identity",
             py::arg("n_jobs")=0,py::arg("validation_ratio")=0.2,py::arg("intercept")=NAN_DOUBLE,
             py::arg("reserved_terms_times_num_x")=100,py::arg("bins")=300,py::arg("verbosity")=0,
             py::arg("max_interaction_level")=1,py::arg("max_interactions")=100000,py::arg("min_observations_in_split")=20,
             py::arg("ineligible_boosting_steps_added")=10,py::arg("max_eligible_terms")=5,
             py::arg("tweedie_power")=1.5,
-            py::arg("validation_tuning_metric")="default"
+            py::arg("validation_tuning_metric")="default",
+            py::arg("quantile")=0.5
             )
         .def("fit", &APLRRegressor::fit,py::arg("X"),py::arg("y"),py::arg("sample_weight")=VectorXd(0),py::arg("X_names")=std::vector<std::string>(),
             py::arg("validation_set_indexes")=std::vector<size_t>(),py::arg("prioritized_predictors_indexes")=std::vector<size_t>(),
@@ -67,24 +69,25 @@ PYBIND11_MODULE(aplr_cpp, m) {
         .def_readwrite("max_training_prediction_or_response",&APLRRegressor::max_training_prediction_or_response)
         .def_readwrite("validation_tuning_metric",&APLRRegressor::validation_tuning_metric)
         .def_readwrite("validation_indexes",&APLRRegressor::validation_indexes)
+        .def_readwrite("quantile",&APLRRegressor::quantile)
         .def(py::pickle(
             [](const APLRRegressor &a) { // __getstate__
                 /* Return a tuple that fully encodes the state of the object */
                 return py::make_tuple(a.m,a.v,a.random_state,a.family,a.n_jobs,a.validation_ratio,a.intercept,a.bins,a.verbosity,
                     a.max_interaction_level,a.max_interactions,a.validation_error_steps,a.term_names,a.term_coefficients,a.terms,a.intercept_steps,
                     a.interactions_eligible,a.min_observations_in_split,a.ineligible_boosting_steps_added,a.max_eligible_terms,
                     a.number_of_base_terms,a.feature_importance,a.link_function,a.tweedie_power,a.min_training_prediction_or_response,a.max_training_prediction_or_response,
-                    a.validation_tuning_metric,a.validation_indexes);
+                    a.validation_tuning_metric,a.validation_indexes,a.quantile);
             },
             [](py::tuple t) { // __setstate__
-                if (t.size() != 28)
+                if (t.size() != 29)
                     throw std::runtime_error("Invalid state!");
 
                 /* Create a new C++ instance */
                 APLRRegressor a(t[0].cast<size_t>(),t[1].cast<double>(),t[2].cast<uint_fast32_t>(),t[3].cast<std::string>(),
                     t[22].cast<std::string>(),t[4].cast<size_t>(),t[5].cast<double>(),
                     t[6].cast<double>(),100,t[7].cast<size_t>(),t[8].cast<size_t>(),t[9].cast<size_t>(),t[10].cast<double>(),t[17].cast<size_t>(),
-                    t[23].cast<double>());
+                    t[23].cast<double>(),t[28].cast<double>());
 
                 a.validation_error_steps=t[11].cast<VectorXd>();
                 a.term_names=t[12].cast<std::vector<std::string>>();

cpp/test ALRRegressor quantile.cpp

@@ -0,0 +1,59 @@
+#include <iostream>
+#include "term.h"
+#include "../dependencies/eigen-master/Eigen/Dense"
+#include <vector>
+#include <numeric>
+#include "APLRRegressor.h"
+#include <cmath>
+
+
+using namespace Eigen;
+
+int main()
+{
+    std::vector<bool> tests;
+    tests.reserve(1000);
+
+    //Model
+    APLRRegressor model{APLRRegressor()};
+    model.m=100;
+    model.v=1.0;
+    model.bins=10;
+    model.n_jobs=1;
+    model.family="quantile";
+    model.verbosity=3;
+    model.max_interaction_level=100;
+    model.max_interactions=30;
+    model.min_observations_in_split=50;
+    model.ineligible_boosting_steps_added=10;
+    model.max_eligible_terms=5;
+    model.quantile=0.5;
+
+    //Data    
+    MatrixXd X_train{load_csv_into_eigen_matrix<MatrixXd>("data/X_train.csv")};
+    MatrixXd X_test{load_csv_into_eigen_matrix<MatrixXd>("data/X_test.csv")}; 
+    VectorXd y_train{load_csv_into_eigen_matrix<MatrixXd>("data/y_train.csv")};    
+    VectorXd y_test{load_csv_into_eigen_matrix<MatrixXd>("data/y_test.csv")}; 
+
+    VectorXd sample_weight{VectorXd::Constant(y_train.size(),1.0)};
+
+    std::cout<<X_train;
+
+    //Fitting
+    //model.fit(X_train,y_train);
+    model.fit(X_train,y_train,sample_weight);
+    //model.fit(X_train,y_train,sample_weight,{},{0,1,2,3,4,5,10,static_cast<size_t>(y_train.size()-1)});
+    std::cout<<"feature importance\n"<<model.feature_importance<<"\n\n";
+
+    VectorXd predictions{model.predict(X_test)};
+    MatrixXd li{model.calculate_local_feature_importance(X_test)};
+
+    //Saving results
+    save_as_csv_file("data/output.csv",predictions);
+
+    std::cout<<predictions.mean()<<"\n\n";
+    tests.push_back(is_approximately_equal(predictions.mean(),23.7883,0.00001));
+
+    //Test summary
+    std::cout<<"\n\nTest summary\n"<<"Passed "<<std::accumulate(tests.begin(),tests.end(),0)<<" out of "<<tests.size()<<" tests.";
+}

examples/train_aplr_cross_validation.py

@@ -31,7 +31,7 @@
 
 #Training model
 param_grid = {"max_interaction_level":[0,1,2,3,100],"min_observations_in_split":[1, 20, 50, 100, 200]}
-family="gaussian" #other available families are binomial, poisson, gamma, tweedie, group_gaussian and mae.
+family="gaussian" #other available families are binomial, poisson, gamma, tweedie, group_gaussian, mae and quantile.
 link_function="identity" #other available link functions are logit and log.
 grid_search_cv = GridSearchCV(APLRRegressor(random_state=random_state,verbosity=1,m=1000,v=0.1,family=family,link_function=link_function),param_grid,cv=5,n_jobs=4,scoring="neg_mean_squared_error")
 grid_search_cv.fit(data_train[predictors].values,data_train[response].values)