From 23651a4a5c20f0af880ca19ad14f72879d1f410c Mon Sep 17 00:00:00 2001 From: Morten Hjorth-Jensen Date: Thu, 29 May 2025 16:53:37 +0200 Subject: [PATCH] Create regression_modelspolybasis.cpp --- .../Regression/regression_modelspolybasis.cpp | 189 ++++++++++++++++++ 1 file changed, 189 insertions(+) create mode 100644 doc/Programs/Regression/regression_modelspolybasis.cpp diff --git a/doc/Programs/Regression/regression_modelspolybasis.cpp b/doc/Programs/Regression/regression_modelspolybasis.cpp new file mode 100644 index 000000000..eff3ed1fd --- /dev/null +++ b/doc/Programs/Regression/regression_modelspolybasis.cpp @@ -0,0 +1,189 @@ +# regression_models_poly_gridsearch.cpp + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; +using namespace Eigen; + +// === Utility Functions === +double mean_squared_error(const VectorXd& y_true, const VectorXd& y_pred) { + return (y_true - y_pred).squaredNorm() / y_true.size(); +} + +double r2_score(const VectorXd& y_true, const VectorXd& y_pred) { + double mean_y = y_true.mean(); + double total = (y_true.array() - mean_y).square().sum(); + double residual = (y_true - y_pred).squaredNorm(); + return 1.0 - residual / total; +} + +void save_csv(const string& filename, const MatrixXd& X, const VectorXd& y_true, const VectorXd& y_pred) { + ofstream file(filename); + file << "X1,X2,...,True Y,Predicted Y\n"; + for (int i = 0; i < X.rows(); ++i) { + for (int j = 0; j < X.cols(); ++j) + file << X(i, j) << (j == X.cols() - 1 ? "," : ","); + file << y_true(i) << "," << y_pred(i) << "\n"; + } + file.close(); +} + +// === Polynomial Expansion === +MatrixXd polynomial_expand(const MatrixXd& X, int degree) { + int n = X.rows(); + int d = X.cols(); + vector terms; + terms.push_back(VectorXd::Ones(n)); // bias + + for (int deg = 1; deg <= degree; ++deg) { + function generate; + generate = [&](int pos, int rem_deg, VectorXi powers) { + if (pos == d) { + if (rem_deg == 0) { + VectorXd term = VectorXd::Ones(n); + for (int i = 0; i < d; ++i) + term = term.array() * X.col(i).array().pow(powers(i)); + terms.push_back(term); + } + return; + } + for (int i = 0; i <= rem_deg; ++i) { + powers(pos) = i; + generate(pos + 1, rem_deg - i, powers); + } + }; + generate(0, deg, VectorXi::Zero(d)); + } + + MatrixXd result(n, terms.size()); + for (int i = 0; i < terms.size(); ++i) + result.col(i) = terms[i]; + return result; +} + +// === Cross-validation === +void cross_validate(const MatrixXd& X, const VectorXd& y, int k, + function fit_func, + function predict_func, + double& avg_mse, double& avg_r2) { + + int n = X.rows(); + vector indices(n); + iota(indices.begin(), indices.end(), 0); + random_shuffle(indices.begin(), indices.end()); + + avg_mse = 0.0; + avg_r2 = 0.0; + + for (int i = 0; i < k; ++i) { + int start = i * n / k; + int end = (i + 1) * n / k; + + vector test_idx(indices.begin() + start, indices.begin() + end); + vector train_idx; + for (int j = 0; j < n; ++j) + if (j < start || j >= end) train_idx.push_back(indices[j]); + + MatrixXd X_train(train_idx.size(), X.cols()); + VectorXd y_train(train_idx.size()); + for (int j = 0; j < train_idx.size(); ++j) { + X_train.row(j) = X.row(train_idx[j]); + y_train(j) = y(train_idx[j]); + } + + MatrixXd X_test(test_idx.size(), X.cols()); + VectorXd y_test(test_idx.size()); + for (int j = 0; j < test_idx.size(); ++j) { + X_test.row(j) = X.row(test_idx[j]); + y_test(j) = y(test_idx[j]); + } + + fit_func(X_train, y_train); + VectorXd y_pred = predict_func(X_test); + avg_mse += mean_squared_error(y_test, y_pred); + avg_r2 += r2_score(y_test, y_pred); + } + + avg_mse /= k; + avg_r2 /= k; +} + +// === Ridge Regression === +class RidgeRegression { +public: + VectorXd weights; + double alpha; + + RidgeRegression(double alpha = 1.0) : alpha(alpha) {} + + void fit(const MatrixXd& X, const VectorXd& y) { + MatrixXd X_bias(X.rows(), X.cols() + 1); + X_bias << MatrixXd::Ones(X.rows(), 1), X; + MatrixXd I = MatrixXd::Identity(X_bias.cols(), X_bias.cols()); + I(0, 0) = 0; + weights = (X_bias.transpose() * X_bias + alpha * I).ldlt().solve(X_bias.transpose() * y); + } + + VectorXd predict(const MatrixXd& X) const { + MatrixXd X_bias(X.rows(), X.cols() + 1); + X_bias << MatrixXd::Ones(X.rows(), 1), X; + return X_bias * weights; + } +}; + +// === Main === +int main() { + // Synthetic data + int n_samples = 150, n_features = 2; + MatrixXd X = MatrixXd::Random(n_samples, n_features); + VectorXd y = 2 + X * VectorXd::LinSpaced(n_features, 1.0, 2.0) + VectorXd::Random(n_samples) * 0.3; + + int best_degree = 1; + double best_alpha = 0.0; + double best_mse = 1e9; + + for (int degree : {1, 2, 3}) { + MatrixXd X_poly = polynomial_expand(X, degree); + for (double alpha : {0.01, 0.1, 1.0, 10.0}) { + RidgeRegression model(alpha); + double avg_mse, avg_r2; + cross_validate(X_poly, y, 5, + [&](const MatrixXd& Xtr, const VectorXd& ytr){ model.fit(Xtr, ytr); }, + [&](const MatrixXd& Xte){ return model.predict(Xte); }, + avg_mse, avg_r2); + + cout << "Degree=" << degree << ", Alpha=" << alpha + << " -> MSE=" << avg_mse << ", R2=" << avg_r2 << endl; + + if (avg_mse < best_mse) { + best_mse = avg_mse; + best_alpha = alpha; + best_degree = degree; + } + } + } + + cout << "\nBest model: Degree=" << best_degree + << ", Alpha=" << best_alpha + << ", MSE=" << best_mse << endl; + + MatrixXd X_poly_best = polynomial_expand(X, best_degree); + RidgeRegression best_model(best_alpha); + best_model.fit(X_poly_best, y); + VectorXd y_pred = best_model.predict(X_poly_best); + + save_csv("predictions_poly_ridge.csv", X, y, y_pred); + + return 0; +} + + +// g++ regression_models_poly_gridsearch.cpp -o poly_model -I /path/to/eigen ./poly_model