10.12.0

Author: mathias-von-ottenbreit

API_REFERENCE_FOR_REGRESSION.md

@@ -50,10 +50,10 @@ Limits 1) the number of terms already in the model that can be considered as int
 Specifies the variance power when ***loss_function*** is "tweedie". Specifies a dispersion parameter when ***loss_function*** is "negative_binomial", "cauchy" or "weibull". 
 
 #### 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" (normalized), "group_mse", "group_mse_by_prediction" and "custom_function". The default is often a choice that fits well with respect to the ***loss_function*** chosen. However, if you want to use ***loss_function*** or ***dispersion_parameter*** as tuning parameters then the default is not suitable. "group_mse" requires that the "group" argument in the ***fit*** method is provided. "group_mse_by_prediction" groups predictions by up to ***group_mse_by_prediction_bins*** groups and calculates groupwise mse. For "custom_function" see ***calculate_custom_validation_error_function*** below.
+Specifies which metric to use for validating the model and tuning ***m***. The model will try to minimize the validation metric. Available options are "default" (using the same methodology as when calculating the training error), "mse", "mae", "negative_gini" (normalized), "group_mse", "group_mse_by_prediction", "neg_top_quantile_mean_response", "bottom_quantile_mean_response" and "custom_function". The default is often a choice that fits well with respect to the ***loss_function*** chosen. However, if you want to use ***loss_function*** or ***dispersion_parameter*** as tuning parameters then the default is not suitable. "group_mse" requires that the "group" argument in the ***fit*** method is provided. "group_mse_by_prediction" groups predictions by up to ***group_mse_by_prediction_bins*** groups and calculates groupwise mse. "neg_top_quantile_mean_response" calculates the negative of the sample weighted mean response for observations with predictions in the top quantile (as specified by the ***quantile*** parameter). For example, if ***quantile*** is 0.95, this metric will be the negative of the sample weighted mean response for the 5% of observations with the highest predictions. "bottom_quantile_mean_response" calculates the sample weighted mean response for observations with predictions in the bottom quantile (as specified by the ***quantile*** parameter). For example, if ***quantile*** is 0.05, this metric will be the sample weighted mean response for the 5% of observations with the lowest predictions. For "custom_function" see ***calculate_custom_validation_error_function*** below.
 
 #### quantile (default = 0.5)
-Specifies the quantile to use when ***loss_function*** is "quantile".
+Specifies the quantile to use when ***loss_function*** is "quantile" or when ***validation_tuning_metric*** is "neg_top_quantile_mean_response" or "bottom_quantile_mean_response".
 
 #### calculate_custom_validation_error_function (default = None)
 A Python function that calculates validation error if ***validation_tuning_metric*** is "custom_function". Example:

cpp/APLRRegressor.h

@@ -143,6 +143,7 @@ class APLRRegressor
     void update_gradient_and_errors();
     void add_new_term(size_t boosting_step);
     void update_coefficient_steps(size_t boosting_step);
+    double calculate_quantile_mean_response(const VectorXd &predictions, bool top_quantile);
     void calculate_and_validate_validation_error(size_t boosting_step);
     double calculate_validation_error(const VectorXd &predictions);
     double calculate_group_mse_by_prediction_validation_error(const VectorXd &predictions);
@@ -190,6 +191,7 @@ class APLRRegressor
     std::string compute_raw_base_term_name(const Term &term, const std::string &X_name);
     void throw_error_if_m_is_invalid();
     bool model_has_not_been_trained();
+    void throw_error_if_quantile_is_invalid();
     std::vector<size_t> compute_relevant_term_indexes(const std::string &unique_term_affiliation);
     std::vector<double> compute_split_points(size_t predictor_index, const std::vector<size_t> &relevant_term_indexes);
     VectorXd compute_contribution_to_linear_predictor_from_specific_terms(const MatrixXd &X, const std::vector<size_t> &term_indexes,
@@ -390,6 +392,7 @@ void APLRRegressor::fit(const MatrixXd &X, const VectorXd &y, const VectorXd &sa
     throw_error_if_loss_function_does_not_exist();
     throw_error_if_link_function_does_not_exist();
     throw_error_if_dispersion_parameter_is_invalid();
+    throw_error_if_quantile_is_invalid();
     throw_error_if_m_is_invalid();
     validate_input_to_fit(X, y, sample_weight, X_names, cv_observations, prioritized_predictors_indexes, monotonic_constraints, group,
                           interaction_constraints, other_data, predictor_learning_rates, predictor_penalties_for_non_linearity,
@@ -606,6 +609,17 @@ void APLRRegressor::throw_error_if_m_is_invalid()
         throw std::runtime_error("The maximum number of boosting steps, m, must be at least 1.");
 }
 
+void APLRRegressor::throw_error_if_quantile_is_invalid()
+{
+    if (loss_function == "quantile" || validation_tuning_metric == "neg_top_quantile_mean_response" || validation_tuning_metric == "bottom_quantile_mean_response")
+    {
+        if (quantile < 0.0 || quantile > 1.0)
+        {
+            throw std::runtime_error("Quantile must be between 0.0 and 1.0.");
+        }
+    }
+}
+
 void APLRRegressor::validate_input_to_fit(const MatrixXd &X, const VectorXd &y, const VectorXd &sample_weight,
                                           const std::vector<std::string> &X_names, const MatrixXi &cv_observations,
                                           const std::vector<size_t> &prioritized_predictors_indexes, const std::vector<int> &monotonic_constraints, const VectorXi &group,
@@ -1720,6 +1734,31 @@ void APLRRegressor::update_coefficient_steps(size_t boosting_step)
     }
 }
 
+double APLRRegressor::calculate_quantile_mean_response(const VectorXd &predictions, bool top_quantile)
+{
+    double quantile_value{calculate_quantile(predictions, quantile, sample_weight_validation)};
+
+    VectorXd predictions_in_quantile;
+    if (top_quantile)
+    {
+        predictions_in_quantile = (predictions.array() >= quantile_value).cast<double>();
+    }
+    else
+    {
+        predictions_in_quantile = (predictions.array() <= quantile_value).cast<double>();
+    }
+
+    VectorXd y_in_quantile{y_validation.array() * predictions_in_quantile.array()};
+    VectorXd weights_in_quantile{sample_weight_validation.array() * predictions_in_quantile.array()};
+
+    double mean_response{calculate_weighted_average(y_in_quantile, weights_in_quantile)};
+
+    if (std::isnan(mean_response))
+        return std::numeric_limits<double>::infinity();
+
+    return mean_response;
+}
+
 void APLRRegressor::calculate_and_validate_validation_error(size_t boosting_step)
 {
     validation_error_steps.col(0)[boosting_step] = calculate_validation_error(predictions_current_validation);
@@ -1784,6 +1823,19 @@ double APLRRegressor::calculate_validation_error(const VectorXd &predictions)
             throw std::runtime_error(error_msg);
         }
     }
+    else if (validation_tuning_metric == "neg_top_quantile_mean_response")
+    {
+        double mean_response{calculate_quantile_mean_response(predictions, true)};
+        if (std::isinf(mean_response))
+        {
+            return mean_response;
+        }
+        return -mean_response;
+    }
+    else if (validation_tuning_metric == "bottom_quantile_mean_response")
+    {
+        return calculate_quantile_mean_response(predictions, false);
+    }
     else
         throw std::runtime_error(validation_tuning_metric + " is an invalid validation_tuning_metric.");
 }

cpp/functions.h

@@ -280,6 +280,27 @@ double calculate_mean_error(const VectorXd &errors, const VectorXd &sample_weigh
     return error;
 }
 
+double calculate_weighted_average(const VectorXd &values, const VectorXd &weights)
+{
+    if (values.size() != weights.size())
+    {
+        throw std::runtime_error("Values and weights must have the same size for weighted average calculation.");
+    }
+    if (values.size() == 0)
+    {
+        return NAN_DOUBLE;
+    }
+
+    double total_weight = weights.sum();
+    if (is_approximately_zero(total_weight))
+    {
+        return NAN_DOUBLE;
+    }
+
+    double weighted_sum = (values.array() * weights.array()).sum();
+    return weighted_sum / total_weight;
+}
+
 double calculate_sum_error(const VectorXd &errors)
 {
     double error{errors.sum()};
@@ -566,4 +587,96 @@ MatrixXd generate_combinations_and_one_additional_column(const std::vector<std::
         }
     }
     return result;
+}
+
+double calculate_quantile(const VectorXd &vector, double quantile, const VectorXd &sample_weight = VectorXd(0))
+{
+    if (quantile < 0.0 || quantile > 1.0)
+    {
+        throw std::runtime_error("Quantile must be between 0.0 and 1.0.");
+    }
+
+    const Eigen::Index n = vector.size();
+    if (n == 0)
+    {
+        return NAN_DOUBLE;
+    }
+
+    VectorXd sample_weight_used;
+    if (sample_weight.size() > 0)
+    {
+        if (sample_weight.size() != n)
+        {
+            throw std::runtime_error("Vector and sample_weight must have the same size.");
+        }
+        sample_weight_used = sample_weight;
+    }
+    else
+    {
+        sample_weight_used = VectorXd::Constant(n, 1.0);
+    }
+
+    if ((sample_weight_used.array() < 0.0).any())
+    {
+        throw std::runtime_error("Sample weights must be non-negative.");
+    }
+
+    double total_weight = sample_weight_used.sum();
+    if (is_approximately_zero(total_weight))
+    {
+        return NAN_DOUBLE;
+    }
+
+    if (n == 1)
+    {
+        return vector[0];
+    }
+
+    std::vector<std::pair<double, double>> weighted_values(n);
+    for (Eigen::Index i = 0; i < n; ++i)
+    {
+        weighted_values[i] = {vector[i], sample_weight_used[i]};
+    }
+
+    std::sort(weighted_values.begin(), weighted_values.end(),
+              [](const auto &a, const auto &b)
+              {
+                  return a.first < b.first;
+              });
+
+    VectorXd quantile_positions(n);
+    double cum_weight = 0.0;
+    for (Eigen::Index i = 0; i < n; ++i)
+    {
+        double current_weight = weighted_values[i].second;
+        cum_weight += current_weight;
+        quantile_positions[i] = (cum_weight - 0.5 * current_weight) / total_weight;
+    }
+
+    auto it = std::upper_bound(quantile_positions.data(), quantile_positions.data() + n, quantile);
+    Eigen::Index upper_index = std::distance(quantile_positions.data(), it);
+
+    if (upper_index == 0)
+    {
+        return weighted_values[0].first;
+    }
+    if (upper_index >= n)
+    {
+        return weighted_values[n - 1].first;
+    }
+
+    Eigen::Index lower_index = upper_index - 1;
+
+    double q_lower = quantile_positions[lower_index];
+    double q_upper = quantile_positions[upper_index];
+    double val_lower = weighted_values[lower_index].first;
+    double val_upper = weighted_values[upper_index].first;
+
+    if (is_approximately_equal(q_lower, q_upper))
+    {
+        return val_lower;
+    }
+
+    double fraction = (quantile - q_lower) / (q_upper - q_lower);
+    return val_lower + fraction * (val_upper - val_lower);
 }

cpp/tests.cpp

@@ -1637,6 +1637,88 @@ class Tests
         tests.push_back(is_approximately_equal(predictions.mean(), 23.646255799722155));
     }
 
+    void test_aplrregressor_neg_top_quantile_mean_response()
+    {
+        // Model
+        APLRRegressor model{APLRRegressor()};
+        model.m = 100;
+        model.v = 1.0;
+        model.bins = 10;
+        model.n_jobs = 1;
+        model.loss_function = "mse";
+        model.validation_tuning_metric = "neg_top_quantile_mean_response";
+        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.8;
+        model.ridge_penalty = 0.0;
+
+        // 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(), 0.5)};
+
+        // Fitting
+        model.fit(X_train, y_train, sample_weight);
+        std::cout << "feature importance\n"
+                  << model.feature_importance << "\n\n";
+
+        VectorXd predictions{model.predict(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.609343969688034));
+    }
+
+    void test_aplrregressor_bottom_quantile_mean_response()
+    {
+        // Model
+        APLRRegressor model{APLRRegressor()};
+        model.m = 100;
+        model.v = 1.0;
+        model.bins = 10;
+        model.n_jobs = 1;
+        model.loss_function = "mse";
+        model.validation_tuning_metric = "bottom_quantile_mean_response";
+        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.2;
+        model.ridge_penalty = 0.0;
+
+        // 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(), 0.5)};
+
+        // Fitting
+        model.fit(X_train, y_train, sample_weight);
+        std::cout << "feature importance\n"
+                  << model.feature_importance << "\n\n";
+
+        VectorXd predictions{model.predict(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.273887245225175));
+    }
+
     void test_aplrregressor_weibull()
     {
         // Model
@@ -2685,6 +2767,8 @@ int main()
     tests.test_aplrregressor_poisson();
     tests.test_aplrregressor_poissongamma();
     tests.test_aplrregressor_quantile();
+    tests.test_aplrregressor_neg_top_quantile_mean_response();
+    tests.test_aplrregressor_bottom_quantile_mean_response();
     tests.test_aplrregressor_weibull();
     tests.test_aplrregressor();
     tests.test_aplr_classifier_multi_class_other_params();

documentation/APLR 10.11.1.pdf renamed to documentation/APLR 10.12.0.pdf

@@ -28,7 +28,7 @@
 
 setuptools.setup(
     name="aplr",
-    version="10.11.1",
+    version="10.12.0",
     description="Automatic Piecewise Linear Regression",
     ext_modules=[sfc_module],
     author="Mathias von Ottenbreit",