Fixed bugs with other model families than gaussian. Added inversegaussian family.

Author: mathias-von-ottenbreit

README.md

@@ -14,7 +14,7 @@ Currently available for Windows and most Linux distributions.
 Please see the two example Python scripts in the examples folder. They cover common use cases, but not all of the functionality in this package. For example, fitting with user-specified observation weights is possible but the example scripts do not use this functionality.
 
 # Sponsorship
-Please consider sponsoring Ottenbreit Data Science by clicking on the Sponsor button. Sufficient funding will enable maintenance of APLR and further development, such as developing a classifier based on APLR.
+Please consider sponsoring Ottenbreit Data Science by clicking on the Sponsor button. Sufficient funding will enable maintenance of APLR and further development.
 
 # API reference
 A thorough API reference will be provided in the future.

cpp/APLRRegressor.h

@@ -79,6 +79,10 @@ class APLRRegressor
     void throw_error_if_family_does_not_exist();
     VectorXd calculate_linear_predictor(const MatrixXd &X);
     void update_linear_predictor_and_predictors();
+    void throw_error_if_response_contains_invalid_values(const VectorXd &y);
+    void throw_error_if_response_is_not_between_0_and_1(const VectorXd &y);
+    void throw_error_if_response_is_negative(const VectorXd &y);
+    void throw_error_if_response_is_not_greater_than_zero(const VectorXd &y);
 
 public:
     //Fields
@@ -191,6 +195,7 @@ void APLRRegressor::validate_input_to_fit(const MatrixXd &X,const VectorXd &y,co
     throw_error_if_matrix_has_nan_or_infinite_elements(y, "y");
     throw_error_if_matrix_has_nan_or_infinite_elements(sample_weight, "sample_weight");
     throw_error_if_validation_set_indexes_has_invalid_indexes(y, validation_set_indexes);
+    throw_error_if_response_contains_invalid_values(y);
 }
 
 void APLRRegressor::throw_error_if_validation_set_indexes_has_invalid_indexes(const VectorXd &y, const std::vector<size_t> &validation_set_indexes)
@@ -205,6 +210,39 @@ void APLRRegressor::throw_error_if_validation_set_indexes_has_invalid_indexes(co
     }
 }
 
+void APLRRegressor::throw_error_if_response_contains_invalid_values(const VectorXd &y)
+{
+    if(family=="logit")
+        throw_error_if_response_is_not_between_0_and_1(y);
+    else if(family=="poisson" || family=="poissongamma")
+        throw_error_if_response_is_negative(y);
+    else if(family=="gamma" || family=="inversegaussian")
+        throw_error_if_response_is_not_greater_than_zero(y);
+}
+
+void APLRRegressor::throw_error_if_response_is_not_between_0_and_1(const VectorXd &y)
+{
+    bool response_is_less_than_zero{(y.array()<0.0).any()};
+    bool response_is_greater_than_one{(y.array()>1.0).any()};
+    if(response_is_less_than_zero || response_is_greater_than_one)
+        throw std::runtime_error("Response values for "+family+" models cannot be less than zero or greater than one.");   
+}
+
+void APLRRegressor::throw_error_if_response_is_negative(const VectorXd &y)
+{
+    bool response_is_less_than_zero{(y.array()<0.0).any()};
+    if(response_is_less_than_zero)
+        throw std::runtime_error("Response values for "+family+" models cannot be less than zero.");   
+}
+
+void APLRRegressor::throw_error_if_response_is_not_greater_than_zero(const VectorXd &y)
+{
+    bool response_is_not_greater_than_zero{(y.array()<=0.0).any()};
+    if(response_is_not_greater_than_zero)
+        throw std::runtime_error("Response values for "+family+" models must be greater than zero.");   
+
+}
+
 void APLRRegressor::define_training_and_validation_sets(const MatrixXd &X,const VectorXd &y,const VectorXd &sample_weight, const std::vector<size_t> &validation_set_indexes)
 {
     //Defining train and validation indexes
@@ -349,6 +387,8 @@ VectorXd APLRRegressor::calculate_neg_gradient_current(const VectorXd &y,const V
         output=(y.array() - predictions_current.array()) / predictions_current.array() / predictions_current.array();
     else if(family=="poissongamma")
         output=(y.array() / predictions_current.array().pow(1.5) - predictions_current.array().pow(-0.5));
+    else if(family=="inversegaussian")
+        output=y.array() / predictions_current.array().pow(3.0) - predictions_current.array().pow(-2.0);
     return output;
 }
 
@@ -648,6 +688,15 @@ void APLRRegressor::select_the_best_term_and_update_errors(size_t boosting_step)
     }
 
     validation_error_steps[boosting_step]=calculate_error(calculate_errors(y_validation,predictions_current_validation,sample_weight_validation,family),sample_weight_validation);
+    bool validation_error_is_invalid{!std::isfinite(validation_error_steps[boosting_step]) || std::isnan(validation_error_steps[boosting_step])};
+    if(validation_error_is_invalid)
+    {
+        abort_boosting=true;
+        std::string warning_message{"Warning: Encountered numerical problems when calculating prediction errors."};
+        if(family=="poisson" || family=="poissongamma" ||family=="gamma" || family=="inversegaussian")
+            warning_message+=" A reason may be too large response values.";
+        std::cout<<warning_message<<"\n";
+    }
 }
 
 void APLRRegressor::update_linear_predictor_and_predictors()
@@ -1012,6 +1061,8 @@ void APLRRegressor::throw_error_if_family_does_not_exist()
         family_exists=true;
     else if(family=="poissongamma")
         family_exists=true;        
+    else if(family=="inversegaussian")
+        family_exists=true;        
     if(!family_exists)
         throw std::runtime_error("Family "+family+" is not available in APLR.");   
 }

cpp/constants.h

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

cpp/functions.h

@@ -71,6 +71,12 @@ VectorXd calculate_poissongamma_errors(const VectorXd &y,const VectorXd &predict
     return errors;
 }
 
+VectorXd calculate_inversegaussian_errors(const VectorXd &y,const VectorXd &predicted)
+{
+    VectorXd errors{y.array().pow(-1.0).array() / 2.0 + y.array()*predicted.array().pow(-2.0) / 2.0 + predicted.array().pow(-1.0) / (-1.0)};
+    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),const std::string &family="gaussian")
 {   
@@ -86,6 +92,8 @@ VectorXd calculate_errors(const VectorXd &y,const VectorXd &predicted,const Vect
         errors=calculate_gamma_errors(y,predicted);
     else if(family=="poissongamma")
         errors=calculate_poissongamma_errors(y,predicted);
+    else if(family=="inversegaussian")
+        errors=calculate_inversegaussian_errors(y,predicted);
     //Adjusting for sample weights if specified
     if(sample_weight.size()>0)
         errors=errors.array()*sample_weight.array();
@@ -127,18 +135,6 @@ 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")
@@ -148,45 +144,11 @@ VectorXd transform_linear_predictor_to_predictions(const VectorXd &linear_predic
         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")
+    else if(family=="poisson" || family=="poissongamma" || family=="gamma" || family=="inversegaussian")
         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/test ALRRegressor gamma.cpp

@@ -51,7 +51,7 @@ int main()
     save_data("data/output.csv",predictions);
 
     std::cout<<predictions.mean()<<"\n\n";
-    tests.push_back(check_if_approximately_equal(predictions.mean(),23.9757,0.00001));
+    tests.push_back(check_if_approximately_equal(predictions.mean(),21.0539,0.00001));
 
     //std::cout<<model.validation_error_steps<<"\n\n";
 

cpp/test ALRRegressor inversegaussian.cpp

@@ -0,0 +1,60 @@
+#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=0.1;
+    model.bins=300;
+    model.n_jobs=0;
+    model.family="inversegaussian";
+    model.verbosity=3;
+    model.max_interaction_level=0;
+    model.max_interactions=1000;
+    model.min_observations_in_split=20;
+    model.ineligible_boosting_steps_added=10;
+    model.max_eligible_terms=5;
+
+    //Data    
+    MatrixXd X_train{load_csv<MatrixXd>("data/X_train.csv")};
+    MatrixXd X_test{load_csv<MatrixXd>("data/X_test.csv")}; 
+    VectorXd y_train{load_csv<MatrixXd>("data/y_train.csv")};    
+    VectorXd y_test{load_csv<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_data("data/output.csv",predictions);
+
+    std::cout<<predictions.mean()<<"\n\n";
+    tests.push_back(check_if_approximately_equal(predictions.mean(),14.8751,0.00001));
+
+    //std::cout<<model.validation_error_steps<<"\n\n";
+
+    //Test summary
+    std::cout<<"\n\nTest summary\n"<<"Passed "<<std::accumulate(tests.begin(),tests.end(),0)<<" out of "<<tests.size()<<" tests.";
+}

cpp/test ALRRegressor poissongamma.cpp

@@ -51,7 +51,7 @@ int main()
     save_data("data/output.csv",predictions);
 
     std::cout<<predictions.mean()<<"\n\n";
-    tests.push_back(check_if_approximately_equal(predictions.mean(),1.89343,0.00001));
+    tests.push_back(check_if_approximately_equal(predictions.mean(),1.88045,0.00001));
 
     //std::cout<<model.validation_error_steps<<"\n\n";
 

examples/train_aplr_cross_validation.py

@@ -28,7 +28,7 @@
 
 #Training model
 param_grid = {"max_interactions":[100000],"max_interaction_level":[0,1,2,3,100],"min_observations_in_split":[1, 20, 50, 100, 200]}
-family="gaussian" #other available families are logit, poisson, gamma and poissongamma
+family="gaussian" #other available families are logit, poisson, poissongamma, gamma and inversegaussian
 grid_search_cv = GridSearchCV(APLRRegressor(random_state=random_state,verbosity=1,m=1000,v=0.1,family=family),param_grid,cv=5,n_jobs=4,scoring="neg_mean_squared_error")
 grid_search_cv.fit(data_train[predictors].values,data_train[response].values)
 best_model:APLRRegressor = grid_search_cv.best_estimator_

examples/train_aplr_validation.py

@@ -32,7 +32,7 @@
 best_validation_result=np.inf
 param_grid=ParameterGrid({"max_interactions":[100000],"max_interaction_level":[0,1,2,3,100],"min_observations_in_split":[1, 20, 50, 100, 200]})
 bestmodel=None
-family="gaussian" #other available families are logit, poisson, gamma and poissongamma
+family="gaussian" #other available families are logit, poisson, poissongamma, gamma and inversegaussian
 for params in param_grid:
     model = APLRRegressor(random_state=random_state,verbosity=3,m=1000,v=0.1,family=family,**params) 
     model.fit(data_train[predictors].values,data_train[response].values,X_names=predictors)

setup.py

@@ -15,7 +15,7 @@
 
 setuptools.setup(
     name='aplr',
-    version='1.1.0',
+    version='1.1.1',
     description='Automatic Piecewise Linear Regression',
     ext_modules=[sfc_module],
     author="Mathias von Ottenbreit",