Множественная регрессия полиномом МНК на C++

Очередной шаблонный класс "Матрица", но достаточно современными средствами C++ и в применении к часто встречающейся задаче построения множественной регрессии, обычно сводящейся к расчёту полинома МНК.

Показаны тесты для МНК "нулевого" (y = const), первого y(x) = a + b * x и второго y(x) = a + b * x + c * x2 порядков, для двух последних тестов данные взяты из английской вики, вычисленные коэффициенты совпали.

Приложение запускалось в консоли актуальной сборки Visual Studio 2019, см. также комментарии в коде.

#include <array>
#include <iostream>

void require(bool condition, const std::string& message) {
 if (condition) return;
 throw std::runtime_error(message);
}

template<typename T, size_t N> //Шаблон для вывода формулы полинома МНК
std::ostream& operator<<(std::ostream& os, const std::array<T, N>& a) {
 auto it = a.cbegin();
 auto end = a.cend();
 os << "y(x) = ";
 if (it != end) {
  os << *it;
  it = std::next(it);
 }
 int pow = 1;
 while (it != end) {
  os << (*it < 0 ? "" : "+") << *it << "*x";
  if (pow > 1) os << "^" << pow;
  pow++;
  it = std::next (it);
 }
 return os;
}

template <size_t RC, size_t CC> //Шаблон для класса "Матрица"
class Matrix {
 std::array<std::array<double, CC>, RC> data;
public:
 Matrix() : data{} {}
 Matrix(std::initializer_list<std::initializer_list<double>> values) {
  size_t rp = 0;
  for (auto row : values) {
   size_t cp = 0;
   for (auto col : row) {
    data[rp][cp] = col;
    cp++;
   }
   rp++;
  }
 }
 double get(size_t row, size_t col) const {
  return data[row][col];
 }
 void set(size_t row, size_t col, double value) {
  data[row][col] = value;
 }
 std::array<double, CC> get(size_t row) {
  return data[row];
 }
 void set(size_t row, const std::array<double, CC>& values) {
  std::copy(values.begin(), values.end(), data[row].begin());
 }

 template <size_t D> //Шаблон для умножения матриц
 Matrix<RC, D> operator* (const Matrix<CC, D>& rhs) const {
  Matrix<RC, D> result;
  for (size_t i = 0; i < RC; i++) {
   for (size_t j = 0; j < D; j++) {
    for (size_t k = 0; k < CC; k++) {
     double prod = get(i, k) * rhs.get(k, j);
     result.set(i, j, result.get(i, j) + prod);
    }
   }
  }
  return result;
 }

 Matrix<CC, RC> transpose() const { //Транспонирование матрицы
  Matrix<CC, RC> trans;
  for (size_t i = 0; i < RC; i++) {
   for (size_t j = 0; j < CC; j++) {
    trans.set(j, i, data[i][j]);
   }
  }
  return trans;
 }

 void toReducedRowEchelonForm() { //Приведение матрицы к ступенчатому виду
  size_t lead = 0;
  for (size_t r = 0; r < RC; r++) {
   if (CC <= lead) return;
   auto i = r;
   while (get(i, lead) == 0.0) {
    i++;
    if (RC == i) {
     i = r;
     lead++;
     if (CC == lead) return;
    }
   }
   auto temp = get(i);
   set(i, get(r));
   set(r, temp);
   if (get(r, lead) != 0.0) {
    auto div = get(r, lead);
    for (size_t j = 0; j < CC; j++) {
     set(r, j, get(r, j) / div);
    }
   }
   for (size_t k = 0; k < RC; k++) {
    if (k != r) {
     auto mult = get(k, lead);
     for (size_t j = 0; j < CC; j++) {
      auto prod = get(r, j) * mult;
      set(k, j, get(k, j) - prod);
     }
    }
   }
   lead++;
  }
 }

 Matrix<RC, RC> inverse() { //Обращение квадратной матрицы
  require(RC == CC, "Not a square matrix");
  Matrix<RC, 2 * RC> aug;
  for (size_t i = 0; i < RC; i++) {
   for (size_t j = 0; j < RC; j++) {
    aug.set(i, j, get(i, j));
   }
   aug.set(i, i + RC, 1.0); //прилепить справа единичную матрицу
  }
  aug.toReducedRowEchelonForm();
  //удалить слева единичную матрицу
  Matrix<RC, RC> inv;
  for (size_t i = 0; i < RC; i++) {
   for (size_t j = RC; j < 2 * RC; j++) {
    inv.set(i, j - RC, aug.get(i, j));
   }
  }
  return inv;
 }
 template <size_t RC, size_t CC>
 friend std::ostream& operator<<(std::ostream&, const Matrix<RC, CC>&);
};

template <size_t RC, size_t CC> //Шаблон для вывода матрицы в поток
std::ostream& operator<<(std::ostream& os, const Matrix<RC, CC>& m) {
 for (size_t i = 0; i < RC; i++) {
  os << '[';
  for (size_t j = 0; j < CC; j++) {
   if (j > 0) {
    os << ", ";
   }
   os << m.get(i, j);
  }
  os << "]" << std::endl;
 }

 return os;
}

template <size_t RC, size_t CC> //Шаблон для основного решения задачи
std::array<double, RC> multiple_regression
(const std::array<double, CC>& y, const Matrix<RC, CC>& x) {
 Matrix<1, CC> tm;
 tm.set(0, y);
 auto cy = tm.transpose();
 auto cx = x.transpose();
 return ((x * cx).inverse() * x * cy).transpose().get(0);
}

//Тесты

void test_power_0() {
 //МНК "нулевого" порядка, то есть y = const
 std::array<double, 5> y{ 1.0, 2.0, 3.0, 4.0, 5.0 };
 Matrix<1, 5> x{ {2.0, 1.0, 3.0, 4.0, 5.0} };
 auto v = multiple_regression(y, x);
 std::cout << v << '\n';
}

void test_power_1() {
 //https://en.wikipedia.org/wiki/Ordinary_least_squares#Solution
 std::array<double, 6> y{ 0.21220, 0.21958, 0.24741, 0.45071, 0.52883, 0.56820 };
 //Матрицу со степенями x для МНК 1 степени зададим статически:
 Matrix<2, 6> x{
  {1.0, 1.0, 1.0, 1.0, 1.0, 1.0 },
  {-0.731354,	-0.707107, -0.615661, 0.052336, 0.309017,	0.438371 }
 };
 auto v = multiple_regression(y, x);
 std::cout << v << '\n';
}

void test_power_2() {
 //https://en.wikipedia.org/wiki/Ordinary_least_squares#Example_with_real_data
 std::array<double, 15> y{ 
  52.21, 53.12, 54.48, 55.84, 57.20, 58.57, 59.93, 61.29, 63.11, 64.47, 
  66.28, 68.10, 69.92, 72.19, 74.46 
 };
 std::array<double, 15> a{ 
  1.47, 1.50, 1.52, 1.55, 1.57, 1.60, 1.63, 1.65, 1.68, 1.70, 
  1.73, 1.75, 1.78, 1.80, 1.83 
 };
 Matrix<3, 15> x;
 //Подготовить матрицу со степенями x для МНК 2 степени
 for (size_t i = 0; i < 15; i++) x.set(0, i, 1.0);
 x.set(1, a);
 for (size_t i = 0; i < 15; i++) x.set(2, i, a[i] * a[i]);
 auto v = multiple_regression(y, x);
 std::cout << v << '\n';
}

int main() {
 test_power_0();
 test_power_1();
 test_power_2();
 return 0;
}

Результаты тестового прогона:

y(x) = 0.981818
y(x) = 0.434783+0.304344*x
y(x) = 128.813-143.162*x+61.9603*x^2

16.02.2022, 12:58 [881 просмотр]