Merge pull request #11 from ottenbreit-data-science/glm

added possibility to train various kinds of glm models

Author: mathias-von-ottenbreit

aplr/aplr.py

@@ -5,11 +5,11 @@
 
 
 class APLRRegressor():
-    def __init__(self, m:int=1000, v:float=0.1, random_state:int=0, loss_function_mse:bool=True, n_jobs:int=0, validation_ratio:float=0.2, intercept:float=np.nan, bins:int=300, max_interaction_level:int=100, max_interactions:int=0, min_observations_in_split:int=20, ineligible_boosting_steps_added:int=10, max_eligible_terms:int=5, verbosity:int=0):
+    def __init__(self, m:int=1000, v:float=0.1, random_state:int=0, family:str="gaussian", n_jobs:int=0, validation_ratio:float=0.2, intercept:float=np.nan, bins:int=300, max_interaction_level:int=100, max_interactions:int=0, min_observations_in_split:int=20, ineligible_boosting_steps_added:int=10, max_eligible_terms:int=5, verbosity:int=0):
         self.m=m
         self.v=v
         self.random_state=random_state
-        self.loss_function_mse=loss_function_mse
+        self.family=family
         self.n_jobs=n_jobs
         self.validation_ratio=validation_ratio
         self.intercept=intercept
@@ -30,7 +30,7 @@ def __set_params_cpp(self):
         self.APLRRegressor.m=self.m
         self.APLRRegressor.v=self.v
         self.APLRRegressor.random_state=self.random_state
-        self.APLRRegressor.loss_function_mse=self.loss_function_mse
+        self.APLRRegressor.family=self.family
         self.APLRRegressor.n_jobs=self.n_jobs
         self.APLRRegressor.validation_ratio=self.validation_ratio
         self.APLRRegressor.intercept=self.intercept
@@ -87,7 +87,7 @@ def get_m(self)->int:
 
     #For sklearn
     def get_params(self, deep=True):
-        return {"m": self.m, "v": self.v,"random_state":self.random_state,"loss_function_mse":self.loss_function_mse,"n_jobs":self.n_jobs,"validation_ratio":self.validation_ratio,"intercept":self.intercept,"bins":self.bins,"max_interaction_level":self.max_interaction_level,"max_interactions":self.max_interactions,"verbosity":self.verbosity,"min_observations_in_split":self.min_observations_in_split,"ineligible_boosting_steps_added":self.ineligible_boosting_steps_added,"max_eligible_terms":self.max_eligible_terms}
+        return {"m": self.m, "v": self.v,"random_state":self.random_state,"family":self.family,"n_jobs":self.n_jobs,"validation_ratio":self.validation_ratio,"intercept":self.intercept,"bins":self.bins,"max_interaction_level":self.max_interaction_level,"max_interactions":self.max_interactions,"verbosity":self.verbosity,"min_observations_in_split":self.min_observations_in_split,"ineligible_boosting_steps_added":self.ineligible_boosting_steps_added,"max_eligible_terms":self.max_eligible_terms}
 
     #For sklearn
     def set_params(self, **parameters):

cpp/APLRRegressor.h

@@ -1,4 +1,5 @@
 #pragma once
 #include <limits>
 
-const double NAN_DOUBLE{ std::numeric_limits<double>::quiet_NaN() };
+const double NAN_DOUBLE{ std::numeric_limits<double>::quiet_NaN() };
+const double SMALL_NEGATIVE_VALUE{-0.001};

cpp/constants.h

@@ -40,38 +40,77 @@ static bool check_if_approximately_zero(TReal a, TReal tolerance = std::numeric_
     return false;
 }
 
+VectorXd calculate_gaussian_errors(const VectorXd &y,const VectorXd &predicted)
+{
+    VectorXd errors{y-predicted};
+    errors=errors.array()*errors.array();
+    return errors;
+}
+
+VectorXd calculate_logit_errors(const VectorXd &y,const VectorXd &predicted)
+{
+    VectorXd errors{-y.array() * predicted.array().log()  -  (1.0-y.array()).array() * (1.0-predicted.array()).log()};
+    return errors;
+}
+
+VectorXd calculate_poisson_errors(const VectorXd &y,const VectorXd &predicted)
+{
+    VectorXd errors{predicted.array() - y.array()*predicted.array().log()};
+    return errors;
+}
+
+VectorXd calculate_gamma_errors(const VectorXd &y,const VectorXd &predicted)
+{
+    VectorXd errors{predicted.array().log() - y.array().log() + y.array()/predicted.array()-1};
+    return errors;
+}
+
+VectorXd calculate_poissongamma_errors(const VectorXd &y,const VectorXd &predicted)
+{
+    VectorXd errors{y.array().pow(0.5).array() / (-0.25) + y.array()*predicted.array().pow(-0.5) / 0.5 + predicted.array().pow(0.5) / 0.5};
+    return errors;
+}
+
 //Computes errors (for each observation) based on error metric for a vector
-VectorXd calculate_errors(const VectorXd &y,const VectorXd &predicted,const VectorXd &sample_weight=VectorXd(0),bool loss_function_mse=true)
+VectorXd calculate_errors(const VectorXd &y,const VectorXd &predicted,const VectorXd &sample_weight=VectorXd(0),const std::string &family="gaussian")
 {   
     //Error per observation before adjustment for sample weights
-    VectorXd residuals{y-predicted};
-    if(loss_function_mse)
-        residuals=residuals.array()*residuals.array();
-    else
-        residuals=residuals.cwiseAbs();
-
+    VectorXd errors;
+    if(family=="gaussian")
+        errors=calculate_gaussian_errors(y,predicted);
+    else if(family=="logit")
+        errors=calculate_logit_errors(y,predicted);
+    else if(family=="poisson")
+        errors=calculate_poisson_errors(y,predicted);
+    else if(family=="gamma")
+        errors=calculate_gamma_errors(y,predicted);
+    else if(family=="poissongamma")
+        errors=calculate_poissongamma_errors(y,predicted);
     //Adjusting for sample weights if specified
     if(sample_weight.size()>0)
-        residuals=residuals.array()*sample_weight.array();
+        errors=errors.array()*sample_weight.array();
 
-    return residuals;
+    return errors;
+}
+
+double calculate_gaussian_error_one_observation(double y,double predicted)
+{
+    double error{y-predicted};
+    error=error*error;
+    return error;
 }
 
 //Computes error for one observation based on error metric
-double calculate_error_one_observation(double y,double predicted,double sample_weight=NAN_DOUBLE,bool loss_function_mse=true)
+double calculate_error_one_observation(double y,double predicted,double sample_weight=NAN_DOUBLE)
 {   
     //Error per observation before adjustment for sample weights
-    double residual{y-predicted};
-    if(loss_function_mse)
-        residual=residual*residual;
-    else
-        residual=abs(residual);
-
+    double error{calculate_gaussian_error_one_observation(y,predicted)};    
+    
     //Adjusting for sample weights if specified
     if(!std::isnan(sample_weight))
-        residual=residual*sample_weight;
-    
-    return residual;
+        error=error*sample_weight;
+
+    return error;
 }
 
 //Computes overall error based on errors from calculate_errors(), returning one value
@@ -88,6 +127,66 @@ double calculate_error(const VectorXd &errors,const VectorXd &sample_weight=Vect
     return error;
 }
 
+VectorXd transform_zero_to_negative(const VectorXd &linear_predictor)
+{
+    VectorXd transformed_linear_predictor{linear_predictor};
+    for (size_t i = 0; i < static_cast<size_t>(transformed_linear_predictor.rows()); ++i)
+    {
+        bool row_is_not_negative{std::isgreaterequal(transformed_linear_predictor[i],0.0)};
+        if(row_is_not_negative)
+            transformed_linear_predictor[i]=SMALL_NEGATIVE_VALUE;
+    }
+    return transformed_linear_predictor;
+}
+
+VectorXd transform_linear_predictor_to_predictions(const VectorXd &linear_predictor, const std::string &family="gaussian")
+{
+    if(family=="gaussian")
+        return linear_predictor;
+    else if(family=="logit")
+    {
+        VectorXd exp_of_linear_predictor{linear_predictor.array().exp()};
+        return exp_of_linear_predictor.array() / (1.0 + exp_of_linear_predictor.array());
+    }
+    else if(family=="poisson")
+        return linear_predictor.array().exp();
+    else if(family=="gamma")
+    {
+        VectorXd transformed_linear_predictor{transform_zero_to_negative(linear_predictor)};
+        return -1/transformed_linear_predictor.array();
+    }
+    else if(family=="poissongamma")
+    {
+        VectorXd transformed_linear_predictor{transform_zero_to_negative(linear_predictor)};
+        return transformed_linear_predictor.array().pow(-2).array() * 4.0;
+    }
+    return VectorXd(0);
+}
+
+double transform_linear_predictor_to_prediction(double linear_predictor, const std::string &family="gaussian")
+{
+    if(family=="gaussian")
+        return linear_predictor;
+    else if(family=="logit")
+    {
+        double exp_of_linear_predictor{std::exp(linear_predictor)};
+        return exp_of_linear_predictor / (1.0 + exp_of_linear_predictor);
+    }
+    else if(family=="poisson")
+        return std::exp(linear_predictor);
+    else if(family=="gamma")
+    {
+        double negative_linear_predictor{std::min(linear_predictor,SMALL_NEGATIVE_VALUE)};
+        return -1/negative_linear_predictor;
+    }
+    else if(family=="poissongamma")
+    {
+        double negative_linear_predictor{std::min(linear_predictor,SMALL_NEGATIVE_VALUE)};
+        return 4.0 * std::pow(negative_linear_predictor,-2);
+    }
+    return NAN_DOUBLE;
+}
+
 //sorts index based on v
 VectorXi sort_indexes_ascending(const VectorXd &v)
 {

cpp/functions.h

@@ -16,7 +16,7 @@ int main()
     model.v=0.5;
     model.bins=300;
     model.n_jobs=0;
-    model.loss_function_mse=true;
+    model.family="gaussian";
     model.verbosity=3;
     model.min_observations_in_split=10;
     //model.max_interaction_level=0;

cpp/main.cpp

@@ -11,8 +11,8 @@ namespace py = pybind11;
 
 PYBIND11_MODULE(aplr_cpp, m) {
     py::class_<APLRRegressor>(m, "APLRRegressor",py::module_local())
-        .def(py::init<int&, double&, int&, bool&,int&,double&,double&,int&,int&,int&,int&,int&,int&,int&,int&>(),
-            py::arg("m")=1000,py::arg("v")=0.1,py::arg("random_state")=0,py::arg("loss_function_mse")=true,
+        .def(py::init<int&, double&, int&, std::string&,int&,double&,double&,int&,int&,int&,int&,int&,int&,int&,int&>(),
+            py::arg("m")=1000,py::arg("v")=0.1,py::arg("random_state")=0,py::arg("family")="gaussian",
             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")=100,py::arg("max_interactions")=0,py::arg("min_observations_in_split")=20,
@@ -40,7 +40,7 @@ PYBIND11_MODULE(aplr_cpp, m) {
         .def_readwrite("max_interactions", &APLRRegressor::max_interactions)
         .def_readwrite("min_observations_in_split", &APLRRegressor::min_observations_in_split)
         .def_readwrite("interactions_eligible", &APLRRegressor::interactions_eligible)
-        .def_readwrite("loss_function_mse", &APLRRegressor::loss_function_mse)
+        .def_readwrite("family", &APLRRegressor::family)
         .def_readwrite("validation_ratio", &APLRRegressor::validation_ratio)
         .def_readwrite("validation_error_steps", &APLRRegressor::validation_error_steps)
         .def_readwrite("n_jobs", &APLRRegressor::n_jobs)
@@ -57,7 +57,7 @@ PYBIND11_MODULE(aplr_cpp, m) {
         .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.loss_function_mse,a.n_jobs,a.validation_ratio,a.intercept,a.bins,a.verbosity,
+                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);
@@ -67,7 +67,7 @@ PYBIND11_MODULE(aplr_cpp, m) {
                     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<bool>(),t[4].cast<size_t>(),t[5].cast<double>(),
+                APLRRegressor a(t[0].cast<size_t>(),t[1].cast<double>(),t[2].cast<uint_fast32_t>(),t[3].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>());
 
                 a.validation_error_steps=t[11].cast<VectorXd>();