67 lines
1.9 KiB
C++
67 lines
1.9 KiB
C++
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
|
// SPDX-License-Identifier: GPL-3.0-only
|
|
|
|
#ifndef JUNGFRAUJOCH_REGRESSION_H
|
|
#define JUNGFRAUJOCH_REGRESSION_H
|
|
|
|
#include <vector>
|
|
#include <utility>
|
|
|
|
#include "../../common/JFJochException.h"
|
|
|
|
struct RegressionResult {
|
|
float intercept;
|
|
float slope;
|
|
float r_square;
|
|
};
|
|
|
|
template <class fptype = double>
|
|
inline RegressionResult regression(std::vector<float> &x, std::vector<float> &y) {
|
|
if (x.size() != y.size())
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
|
"Regression: mismatch between input parameter size");
|
|
|
|
if (x.size() < 3)
|
|
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
|
"Regression: fitting less than 3 points is non-sense");
|
|
|
|
fptype x_sum = 0.0;
|
|
fptype y_sum = 0.0;
|
|
fptype xy_sum = 0.0;
|
|
fptype x2_sum = 0.0;
|
|
|
|
for (int i = 0; i < x.size(); i++) {
|
|
x_sum += x[i];
|
|
y_sum += y[i];
|
|
x2_sum += x[i] * x[i];
|
|
xy_sum += x[i] * y[i];
|
|
}
|
|
|
|
fptype denominator = x.size() * x2_sum - x_sum * x_sum;
|
|
fptype intercept = (y_sum * x2_sum - x_sum * xy_sum) / denominator;
|
|
fptype slope = (x.size() * xy_sum - x_sum * y_sum) / denominator;
|
|
|
|
fptype y_mean = y_sum / y.size();
|
|
fptype ss_tot = 0.0;
|
|
fptype ss_reg = 0.0;
|
|
|
|
for (int i = 0; i < y.size(); i++) {
|
|
fptype pred = slope * x[i] + intercept;
|
|
ss_tot += (y[i] - y_mean) * (y[i] - y_mean);
|
|
ss_reg += (pred - y_mean) * (pred - y_mean);
|
|
}
|
|
|
|
fptype r_square = 0.0;
|
|
|
|
if (ss_tot != 0.0)
|
|
r_square = ss_reg / ss_tot;
|
|
|
|
return {
|
|
.intercept = static_cast<float>(intercept),
|
|
.slope = static_cast<float>(slope),
|
|
.r_square = static_cast<float>(r_square)
|
|
};
|
|
};
|
|
|
|
#endif //JUNGFRAUJOCH_REGRESSION_H
|