intercept methodology

Author: mathias-von-ottenbreit

API_REFERENCE_FOR_REGRESSION.md

@@ -181,11 +181,6 @@ The index of the term selected. So ***0*** is the first term, ***1*** is the sec
 ***Returns the regression coefficient of the intercept term.***
 
 
-## Method: get_intercept_steps()
-
-***Returns a numpy vector containing the regression coefficients of the intercept term by boosting step.***
-
-
 ## Method: get_optimal_m()
 
 ***Returns the number of boosting steps in the model (the value that minimized validation error).***

aplr/aplr.py

@@ -87,9 +87,6 @@ def get_feature_importance(self)->npt.ArrayLike:
     def get_intercept(self)->float:
         return self.APLRRegressor.get_intercept()
 
-    def get_intercept_steps(self)->npt.ArrayLike:
-        return self.APLRRegressor.get_intercept_steps()
-
     def get_optimal_m(self)->int:
         return self.APLRRegressor.get_optimal_m()
 

cpp/APLRRegressor.h

@@ -50,7 +50,6 @@ class APLRRegressor
     std::set<int> unique_groups_train;
     std::set<int> unique_groups_validation;
     std::vector<int> interaction_constraints;
-    bool update_intercept_only_once;
 
     void validate_input_to_fit(const MatrixXd &X,const VectorXd &y,const VectorXd &sample_weight,const std::vector<std::string> &X_names, 
         const std::vector<size_t> &validation_set_indexes, const std::vector<size_t> &prioritized_predictors_indexes,
@@ -78,7 +77,7 @@ class APLRRegressor
     void add_necessary_given_terms_to_interaction(Term &interaction, Term &existing_model_term);
     void find_sorted_indexes_for_errors_for_interactions_to_consider();
     void add_promising_interactions_and_select_the_best_one();
-    void update_intercept(size_t boosting_step);
+    void update_intercept();
     void select_the_best_term_and_update_errors(size_t boosting_step, bool not_evaluating_prioritized_predictors=true);
     void update_terms(size_t boosting_step);
     void update_gradient_and_errors();
@@ -128,7 +127,6 @@ class APLRRegressor
     std::vector<std::string> term_names;
     VectorXd term_coefficients;
     size_t max_interaction_level;
-    VectorXd intercept_steps;
     size_t max_interactions;
     size_t interactions_eligible;
     VectorXd validation_error_steps;
@@ -165,7 +163,6 @@ class APLRRegressor
     VectorXd get_validation_error_steps();
     VectorXd get_feature_importance();
     double get_intercept();
-    VectorXd get_intercept_steps();
     size_t get_optimal_m();
     std::string get_validation_tuning_metric();
     std::vector<size_t> get_validation_indexes();
@@ -178,7 +175,7 @@ APLRRegressor::APLRRegressor(size_t m,double v,uint_fast32_t random_state,std::s
         reserved_terms_times_num_x{reserved_terms_times_num_x},intercept{intercept},m{m},v{v},
         loss_function{loss_function},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)},
+        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},dispersion_parameter{dispersion_parameter},min_training_prediction_or_response{NAN_DOUBLE},
         max_training_prediction_or_response{NAN_DOUBLE}, validation_tuning_metric{validation_tuning_metric},
@@ -191,8 +188,8 @@ APLRRegressor::APLRRegressor(const APLRRegressor &other):
     loss_function{other.loss_function},link_function{other.link_function},validation_ratio{other.validation_ratio},
     n_jobs{other.n_jobs},random_state{other.random_state},bins{other.bins},
     verbosity{other.verbosity},term_names{other.term_names},term_coefficients{other.term_coefficients},
-    max_interaction_level{other.max_interaction_level},intercept_steps{other.intercept_steps},
-    max_interactions{other.max_interactions},interactions_eligible{other.interactions_eligible},validation_error_steps{other.validation_error_steps},
+    max_interaction_level{other.max_interaction_level},max_interactions{other.max_interactions},
+    interactions_eligible{other.interactions_eligible},validation_error_steps{other.validation_error_steps},
     min_observations_in_split{other.min_observations_in_split},ineligible_boosting_steps_added{other.ineligible_boosting_steps_added},
     max_eligible_terms{other.max_eligible_terms},number_of_base_terms{other.number_of_base_terms},
     feature_importance{other.feature_importance},dispersion_parameter{other.dispersion_parameter},min_training_prediction_or_response{other.min_training_prediction_or_response},
@@ -289,6 +286,12 @@ void APLRRegressor::throw_error_if_dispersion_parameter_is_invalid()
     }
 }
 
+void APLRRegressor::throw_error_if_m_is_invalid()
+{
+    if(m<1)
+        throw std::runtime_error("The maximum number of boosting steps, m, must be at least 1.");
+}
+
 void APLRRegressor::validate_input_to_fit(const MatrixXd &X,const VectorXd &y,const VectorXd &sample_weight,
     const std::vector<std::string> &X_names, const std::vector<size_t> &validation_set_indexes, 
     const std::vector<size_t> &prioritized_predictors_indexes, const std::vector<int> &monotonic_constraints, const VectorXi &group,
@@ -546,21 +549,7 @@ void APLRRegressor::initialize(const std::vector<size_t> &prioritized_predictors
     terms.clear();
     terms.reserve(X_train.cols()*reserved_terms_times_num_x);
 
-    if(loss_function == "group_mse")
-    {
-        update_intercept_only_once = true;
-        if(sample_weight_train.size()==0)
-            intercept = y_train.mean();
-        else
-            intercept = (y_train.array()*sample_weight_train.array()).sum()/sample_weight_train.array().sum();
-    }
-    else
-    {
-        update_intercept_only_once = false;
-        intercept=0;
-    }
-    intercept_steps=VectorXd::Constant(m, intercept);
-
+    double initial_prediction{0.0};
 
     terms_eligible_current.reserve(X_train.cols()*reserved_terms_times_num_x);
     size_t X_train_cols{static_cast<size_t>(X_train.cols())};
@@ -602,9 +591,9 @@ void APLRRegressor::initialize(const std::vector<size_t> &prioritized_predictors
         }
     }
 
-    linear_predictor_current=VectorXd::Constant(y_train.size(),intercept);
+    linear_predictor_current=VectorXd::Constant(y_train.size(),initial_prediction);
     linear_predictor_null_model=linear_predictor_current;
-    linear_predictor_current_validation=VectorXd::Constant(y_validation.size(),intercept);
+    linear_predictor_current_validation=VectorXd::Constant(y_validation.size(),initial_prediction);
     predictions_current=transform_linear_predictor_to_predictions(linear_predictor_current,link_function);
     predictions_current_validation=transform_linear_predictor_to_predictions(linear_predictor_current_validation,link_function);
 
@@ -717,17 +706,31 @@ VectorXd APLRRegressor::differentiate_predictions()
 void APLRRegressor::execute_boosting_steps()
 {
     abort_boosting = false;
-    for (size_t boosting_step = 0; boosting_step < m; ++boosting_step)
+    update_intercept();
+    for (size_t boosting_step = 1; boosting_step < m; ++boosting_step)
     {
         execute_boosting_step(boosting_step);
         if(abort_boosting) break;
     }
 }
 
+void APLRRegressor::update_intercept()
+{
+    double intercept_update;
+    if(sample_weight_train.size()==0)
+        intercept=neg_gradient_current.mean();
+    else
+        intercept=(neg_gradient_current.array()*sample_weight_train.array()).sum()/sample_weight_train.array().sum();
+    linear_predictor_update=VectorXd::Constant(neg_gradient_current.size(),intercept);
+    linear_predictor_update_validation=VectorXd::Constant(y_validation.size(),intercept);
+    update_linear_predictor_and_predictions();
+    update_gradient_and_errors();
+    calculate_and_validate_validation_error(0);
+    print_summary_after_boosting_step(0);
+}
+
 void APLRRegressor::execute_boosting_step(size_t boosting_step)
 {
-    if(!update_intercept_only_once)
-        update_intercept(boosting_step);
     bool prioritize_predictors{!abort_boosting && prioritized_predictors_indexes.size()>0};
     if(prioritize_predictors)
     {
@@ -758,25 +761,6 @@ void APLRRegressor::execute_boosting_step(size_t boosting_step)
     print_summary_after_boosting_step(boosting_step);
 }
 
-void APLRRegressor::update_intercept(size_t boosting_step)
-{
-    double intercept_update;
-    if(sample_weight_train.size()==0)
-        intercept_update=v*neg_gradient_current.mean();
-    else
-        intercept_update=v*(neg_gradient_current.array()*sample_weight_train.array()).sum()/sample_weight_train.array().sum();
-    linear_predictor_update=VectorXd::Constant(neg_gradient_current.size(),intercept_update);
-    linear_predictor_update_validation=VectorXd::Constant(y_validation.size(),intercept_update);
-    update_linear_predictor_and_predictions();
-    update_gradient_and_errors();
-    calculate_and_validate_validation_error(boosting_step);
-    if(!abort_boosting)
-    {
-        intercept+=intercept_update;
-        intercept_steps[boosting_step]=intercept;
-    }
-}
-
 void APLRRegressor::update_linear_predictor_and_predictions()
 {
     linear_predictor_current+=linear_predictor_update;
@@ -1219,13 +1203,6 @@ void APLRRegressor::print_summary_after_boosting_step(size_t boosting_step)
 
 void APLRRegressor::update_coefficients_for_all_steps()
 {
-    for (size_t j = 0; j < m; ++j)
-    {
-        bool fill_down_coefficient_steps{j>0 && is_approximately_zero(intercept_steps[j]) && !is_approximately_zero(intercept_steps[j-1])};
-        if(fill_down_coefficient_steps)
-            intercept_steps[j]=intercept_steps[j-1];
-    }
-
     for (size_t i = 0; i < terms.size(); ++i)
     {
         for (size_t j = 0; j < m; ++j)
@@ -1249,7 +1226,6 @@ void APLRRegressor::find_optimal_m_and_update_model_accordingly()
 {
     size_t best_boosting_step_index;
     validation_error_steps.minCoeff(&best_boosting_step_index);
-    intercept=intercept_steps[best_boosting_step_index];
     for (size_t i = 0; i < terms.size(); ++i)
     {
         terms[i].coefficient = terms[i].coefficient_steps[best_boosting_step_index];
@@ -1274,10 +1250,6 @@ void APLRRegressor::revert_scaling_if_using_log_link_function()
     {
         y_train/=scaling_factor_for_log_link_function;
         intercept+=std::log(1/scaling_factor_for_log_link_function);
-        for (Eigen::Index i = 0; i < intercept_steps.size(); ++i)
-        {
-            intercept_steps[i]+=std::log(1/scaling_factor_for_log_link_function);
-        }
     }
 }
 
@@ -1451,13 +1423,13 @@ VectorXd APLRRegressor::predict(const MatrixXd &X, bool cap_predictions_to_minma
 
 VectorXd APLRRegressor::calculate_linear_predictor(const MatrixXd &X)
 {
-    VectorXd predictions{VectorXd::Constant(X.rows(),intercept)};
+    VectorXd linear_predictor{VectorXd::Constant(X.rows(),intercept)};
     for (size_t i = 0; i < terms.size(); ++i)
     {
         VectorXd contrib{terms[i].calculate_contribution_to_linear_predictor(X)};
-        predictions+=contrib;
+        linear_predictor+=contrib;
     }
-    return predictions;    
+    return linear_predictor;    
 }
 
 void APLRRegressor::cap_predictions_to_minmax_in_training(VectorXd &predictions)
@@ -1531,11 +1503,6 @@ double APLRRegressor::get_intercept()
     return intercept;
 }
 
-VectorXd APLRRegressor::get_intercept_steps()
-{
-    return intercept_steps;
-}
-
 size_t APLRRegressor::get_optimal_m()
 {
     return m_optimal;
@@ -1549,10 +1516,4 @@ std::string APLRRegressor::get_validation_tuning_metric()
 std::vector<size_t> APLRRegressor::get_validation_indexes()
 {
     return validation_indexes;
-}
-
-void APLRRegressor::throw_error_if_m_is_invalid()
-{
-    if(m<1)
-        throw std::runtime_error("The maximum number of boosting steps, m, must be at least 1.");   
 }

cpp/pythonbinding.cpp

@@ -39,12 +39,10 @@ PYBIND11_MODULE(aplr_cpp, m) {
         .def("get_validation_error_steps", &APLRRegressor::get_validation_error_steps)
         .def("get_feature_importance", &APLRRegressor::get_feature_importance)
         .def("get_intercept", &APLRRegressor::get_intercept)
-        .def("get_intercept_steps", &APLRRegressor::get_intercept_steps)
         .def("get_optimal_m", &APLRRegressor::get_optimal_m)
         .def("get_validation_tuning_metric", &APLRRegressor::get_validation_tuning_metric)
         .def("get_validation_indexes", &APLRRegressor::get_validation_indexes)
         .def_readwrite("intercept", &APLRRegressor::intercept)
-        .def_readwrite("intercept_steps", &APLRRegressor::intercept_steps)
         .def_readwrite("m", &APLRRegressor::m)
         .def_readwrite("m_optimal", &APLRRegressor::m_optimal)
         .def_readwrite("v", &APLRRegressor::v)
@@ -77,36 +75,35 @@ PYBIND11_MODULE(aplr_cpp, m) {
             [](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,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.max_interaction_level,a.max_interactions,a.validation_error_steps,a.term_names,a.term_coefficients,a.terms,
                     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.dispersion_parameter,a.min_training_prediction_or_response,a.max_training_prediction_or_response,
                     a.validation_tuning_metric,a.validation_indexes,a.quantile,a.m_optimal);
             },
             [](py::tuple t) { // __setstate__
-                if (t.size() != 30)
+                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[28].cast<double>());
+                    t[21].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[16].cast<size_t>(),
+                    t[22].cast<double>(),t[27].cast<double>());
 
                 a.validation_error_steps=t[11].cast<VectorXd>();
                 a.term_names=t[12].cast<std::vector<std::string>>();
                 a.term_coefficients=t[13].cast<VectorXd>();
                 a.terms=t[14].cast<std::vector<Term>>();
-                a.intercept_steps=t[15].cast<VectorXd>();
-                a.interactions_eligible=t[16].cast<size_t>();
-                a.ineligible_boosting_steps_added=t[18].cast<size_t>();
-                a.max_eligible_terms=t[19].cast<size_t>();
-                a.number_of_base_terms=t[20].cast<size_t>();
-                a.feature_importance=t[21].cast<VectorXd>();
-                a.min_training_prediction_or_response=t[24].cast<double>();
-                a.max_training_prediction_or_response=t[25].cast<double>();
-                a.validation_tuning_metric=t[26].cast<std::string>();
-                a.validation_indexes=t[27].cast<std::vector<size_t>>();
-                a.m_optimal=t[29].cast<size_t>();
+                a.interactions_eligible=t[15].cast<size_t>();
+                a.ineligible_boosting_steps_added=t[17].cast<size_t>();
+                a.max_eligible_terms=t[18].cast<size_t>();
+                a.number_of_base_terms=t[19].cast<size_t>();
+                a.feature_importance=t[20].cast<VectorXd>();
+                a.min_training_prediction_or_response=t[23].cast<double>();
+                a.max_training_prediction_or_response=t[24].cast<double>();
+                a.validation_tuning_metric=t[25].cast<std::string>();
+                a.validation_indexes=t[26].cast<std::vector<size_t>>();
+                a.m_optimal=t[28].cast<size_t>();
 
                 return a;
             }

cpp/test ALRRegressor cauchy group_mse validation.cpp

@@ -55,7 +55,7 @@ int main()
     save_as_csv_file("data/output.csv",predictions);
 
     std::cout<<predictions.mean()<<"\n\n";
-    tests.push_back(is_approximately_equal(predictions.mean(),20.2462,0.00001));
+    tests.push_back(is_approximately_equal(predictions.mean(),12.764,0.00001));
 
     //Test summary
     std::cout<<"\n\nTest summary\n"<<"Passed "<<std::accumulate(tests.begin(),tests.end(),0)<<" out of "<<tests.size()<<" tests.";

cpp/test ALRRegressor cauchy.cpp

@@ -52,7 +52,7 @@ int main()
     save_as_csv_file("data/output.csv",predictions);
 
     std::cout<<predictions.mean()<<"\n\n";
-    tests.push_back(is_approximately_equal(predictions.mean(),19.7161,0.00001));
+    tests.push_back(is_approximately_equal(predictions.mean(),12.108,0.00001));
 
     //Test summary
     std::cout<<"\n\nTest summary\n"<<"Passed "<<std::accumulate(tests.begin(),tests.end(),0)<<" out of "<<tests.size()<<" tests.";

cpp/test ALRRegressor gamma rank unweighted.cpp

@@ -53,7 +53,7 @@ int main()
     save_as_csv_file("data/output.csv",predictions);
 
     std::cout<<predictions.mean()<<"\n\n";
-    tests.push_back(is_approximately_equal(predictions.mean(),23.6581,0.00001));
+    tests.push_back(is_approximately_equal(predictions.mean(),23.6509,0.00001));
 
     //Test summary
     std::cout<<"\n\nTest summary\n"<<"Passed "<<std::accumulate(tests.begin(),tests.end(),0)<<" out of "<<tests.size()<<" tests.";

cpp/test ALRRegressor gamma rank.cpp

@@ -53,7 +53,7 @@ int main()
     save_as_csv_file("data/output.csv",predictions);
 
     std::cout<<predictions.mean()<<"\n\n";
-    tests.push_back(is_approximately_equal(predictions.mean(),23.6581,0.00001));
+    tests.push_back(is_approximately_equal(predictions.mean(),23.6509,0.00001));
 
     //Test summary
     std::cout<<"\n\nTest summary\n"<<"Passed "<<std::accumulate(tests.begin(),tests.end(),0)<<" out of "<<tests.size()<<" tests.";

cpp/test ALRRegressor gamma.cpp

@@ -53,7 +53,7 @@ int main()
     save_as_csv_file("data/output.csv",predictions);
 
     std::cout<<predictions.mean()<<"\n\n";
-    tests.push_back(is_approximately_equal(predictions.mean(),23.6503,0.00001));
+    tests.push_back(is_approximately_equal(predictions.mean(),23.6458,0.00001));
 
     //Test summary
     std::cout<<"\n\nTest summary\n"<<"Passed "<<std::accumulate(tests.begin(),tests.end(),0)<<" out of "<<tests.size()<<" tests.";

cpp/test ALRRegressor int constr.cpp

@@ -52,7 +52,7 @@ int main()
     save_as_csv_file("data/output.csv",predictions);
 
     std::cout<<predictions.mean()<<"\n\n";
-    tests.push_back(is_approximately_equal(predictions.mean(),23.4597,0.00001));
+    tests.push_back(is_approximately_equal(predictions.mean(),23.5213,0.00001));
 
     //std::cout<<model.validation_error_steps<<"\n\n";