C++ source code to calculate a simple regression curve(5d).

calc.cpp

#include "calc.hpp"

/*
 * @brief       単回帰曲線（5次）の計算
 *
 * @param[ref]  a (double)
 * @param[ref]  b (double)
 * @param[ref]  c (double)
 * @param[ref]  d (double)
 * @param[ref]  e (double)
 * @param[ref]  f (double)
 * @return      真偽(bool)
 * @retval      true  成功
 * @retval      false 失敗
 */
bool Calc::reg_curve_5d(
    double& a, double& b, double& c, double& d, double& e, double& f) {
  unsigned int i;      // loop インデックス
  double s_x   = 0.0;  // sum(x)
  double s_x2  = 0.0;  // sum(xx)
  double s_x3  = 0.0;  // sum(x^3)
  double s_x4  = 0.0;  // sum(x^4)
  double s_x5  = 0.0;  // sum(x^5)
  double s_x6  = 0.0;  // sum(x^6)
  double s_x7  = 0.0;  // sum(x^7)
  double s_x8  = 0.0;  // sum(x^8)
  double s_x9  = 0.0;  // sum(x^9)
  double s_x10 = 0.0;  // sum(x^10)
  double s_y   = 0.0;  // sum(y)
  double s_xy  = 0.0;  // sum(xy)
  double s_x2y = 0.0;  // sum(xxy)
  double s_x3y = 0.0;  // sum(x^3y)
  double s_x4y = 0.0;  // sum(x^4y)
  double s_x5y = 0.0;  // sum(x^5y)
  double x     = 0.0;  // x    計算用
  double x2    = 0.0;  // xx   計算用
  double x3    = 0.0;  // x^3  計算用
  double x4    = 0.0;  // x^4  計算用
  double x5    = 0.0;  // x^5  計算用
  double x6    = 0.0;  // x^6  計算用
  double x7    = 0.0;  // x^7  計算用
  double x8    = 0.0;  // x^8  計算用
  double x9    = 0.0;  // x^9  計算用
  double x10   = 0.0;  // x^10 計算用
  double y     = 0.0;  // y    計算用

  try {
    cnt = data.size();
    for (i = 0; i < cnt; i++) {
      x  = data[i][0];
      y  = data[i][1];
      x2  = x  * x;
      x3  = x2 * x;
      x4  = x3 * x;
      x5  = x4 * x;
      x6  = x5 * x;
      x7  = x6 * x;
      x8  = x7 * x;
      x9  = x8 * x;
      x10 = x9 * x;
      s_x   += x;
      s_x2  += x2;
      s_x3  += x3;
      s_x4  += x4;
      s_x5  += x5;
      s_x6  += x6;
      s_x7  += x7;
      s_x8  += x8;
      s_x9  += x9;
      s_x10 += x10;
      s_y   += y;
      s_xy  += x  * y;
      s_x2y += x2 * y;
      s_x3y += x3 * y;
      s_x4y += x4 * y;
      s_x5y += x5 * y;
    }

    // 行列1行目
    mtx.push_back({(double)cnt, s_x, s_x2, s_x3, s_x4, s_x5, s_y});
    // 行列2行目
    mtx.push_back({s_x, s_x2, s_x3, s_x4, s_x5, s_x6, s_xy});
    // 行列3行目
    mtx.push_back({s_x2, s_x3, s_x4, s_x5, s_x6, s_x7, s_x2y});
    // 行列4行目
    mtx.push_back({s_x3, s_x4, s_x5, s_x6, s_x7, s_x8, s_x3y});
    // 行列5行目
    mtx.push_back({s_x4, s_x5, s_x6, s_x7, s_x8, s_x9, s_x4y});
    // 行列6行目
    mtx.push_back({s_x5, s_x6, s_x7, s_x8, s_x9, s_x10, s_x5y});

    // 計算（ガウスの消去法）
    if (!solve_ge(mtx)) {
      std::cout << "[ERROR] Failed to solve by the Gauss-Ellimination method!"
                << std::endl;
      return false;
    }

    // a, ..., f
    a = mtx[0][6];
    b = mtx[1][6];
    c = mtx[2][6];
    d = mtx[3][6];
    e = mtx[4][6];
    f = mtx[5][6];
  } catch (...) {
    return false;  // 計算失敗
  }

  return true;  // 計算成功
}

/*
 * @brief      ピボット選択
 *
 * @param[ref] 行列（配列） mtx (vector<vector<double>>)
 * @param[in]  対象行
 * @return     true|false
 */
bool Calc::pivot(std::vector<std::vector<double>>& mtx, unsigned int k) {
  unsigned int n;      // 行数
  unsigned int pv;     // ピボット行
  unsigned int i;      // loop index
  double       v_max;  // k 列絶対値の最大値

  try {
    n = mtx.size();
    pv = k;
    v_max = std::fabs(mtx[k][k]);
    for (i = k + 1; i < n; ++i) {
      if (std::fabs(mtx[i][k]) > v_max) {
        pv = i;
        v_max = std::fabs(mtx[i][k]);
      }
    }
    if (k != pv) { swap(mtx[k], mtx[pv]); }
  } catch (...) {
    return false;
  }

  return true;
}

/**
 * @brief       ガウスの消去法（ピボット選択）
 *
 * @param[ref]  行列（配列） mtx (double)
 * @return      真偽(bool)
 * @retval      true  成功
 * @retval      false 失敗
 */
bool Calc::solve_ge(std::vector<std::vector<double>>& mtx) {
  int i;     // loop インデックス
  int j;     // loop インデックス
  int k;     // loop インデックス
  int n;     // 元（行）の数
  double d;  // 計算用

  try {
    n = (int)mtx.size();

    // 前進消去
    for (k = 0; k < n - 1; k++) {
      if (!pivot(mtx, k)) {
        std::cout << "[ERROR] Failed to pivot!" << std::endl;
        return false;
      }
      for (i = k + 1; i < n; i++) {
        d = mtx[i][k] / mtx[k][k];
        for (j = k + 1; j <= n; j++)
          mtx[i][j] -= mtx[k][j] * d;
      }
    }

    // 後退代入
    for (i = n - 1; i >= 0; i--) {
      d = mtx[i][n];
      for (j = i + 1; j < n; j++)
        d -= mtx[i][j] * mtx[j][n];
      mtx[i][n] = d / mtx[i][i];
    }
  } catch (...) {
    return false;  // 計算失敗
  }

  return true;  // 計算成功
}

calc.hpp

#ifndef REGRESSION_CURVE_5D_CALC_HPP_
#define REGRESSION_CURVE_5D_CALC_HPP_

#include <cmath>
#include <iostream>
#include <sstream>
#include <vector>

class Calc {
  std::vector<std::vector<double>> data;             // 元データ
  std::vector<std::vector<double>> mtx;              // 計算用行列
  bool pivot(std::vector<std::vector<double>>&, unsigned int);
                                                     // ピボット選択
  bool solve_ge(std::vector<std::vector<double>>&);  // ガウスの消去法

public:
  Calc(std::vector<std::vector<double>>& data) : data(data) {}
  unsigned int cnt;                                  // データ件数
  bool reg_curve_5d(double&, double&, double&, double&, double&, double&);
                                                     // 単回帰曲線（5次）の計算
};

#endif

file.cpp

#include "file.hpp"

#include <iostream>
#include <sstream>
#include <string>
#include <vector>

bool File::get_text(std::vector<std::vector<double>>& data) {
  try {
    // ファイル OPEN
    std::ifstream ifs(f_data);
    if (!ifs.is_open()) return false;  // 読み込み失敗

    // ファイル READ
    std::string buf;                   // 1行分バッファ
    while (getline(ifs, buf)) {
      std::vector<double> rec;         // 1行分ベクタ
      std::istringstream iss(buf);     // 文字列ストリーム
      std::string field;               // 1列分文字列

      // 1行分文字列を1行分ベクタに追加
      double x, y;
      while (iss >> x >> y) {
        rec.push_back(x);
        rec.push_back(y);
      }

      // １行分ベクタを data ベクタに追加
      if (rec.size() != 0) data.push_back(rec);
    }
  } catch (...) {
      std::cerr << "EXCEPTION!" << std::endl;
      return false;
  }

  return true;  // 読み込み成功
}

file.hpp

#ifndef REGRESSION_CURVE_5D_FILE_HPP_
#define REGRESSION_CURVE_5D_FILE_HPP_

#include <fstream>
#include <string>
#include <vector>

class File {
  std::string f_data;

public:
  File(std::string f_data) : f_data(f_data) {}
  bool get_text(std::vector<std::vector<double>>&);
};

#endif

regression_curve_5d.cpp

/***********************************************************
  単回帰曲線（5次回帰モデル）計算
  : y = a + b * x + c * x^2 + d * x^3 + e * x^4 + f * x^5
  : 連立方程式をガウス（ピボット選択）の消去法で解く方法

    DATE          AUTHOR          VERSION
    2022.06.20    mk-mode.com     1.00 新規作成

  Copyright(C) 2022 mk-mode.com All Rights Reserved.
***********************************************************/
#include "calc.hpp"
#include "file.hpp"

#include <cstdlib>   // for EXIT_XXXX
#include <iomanip>   // for setprecision
#include <iostream>
#include <string>
#include <vector>

int main(int argc, char* argv[]) {
  std::string f_data;                     // データファイル名
  std::vector<std::vector<double>> data;  // データ配列
  std::size_t i;                          // loop インデックス
  double a;                               // 係数 a
  double b;                               // 係数 b
  double c;                               // 係数 c
  double d;                               // 係数 d
  double e;                               // 定数 e
  double f;                               // 定数 f

  try {
    // コマンドライン引数のチェック
    if (argc != 2) {
      std::cerr << "[ERROR] Number of arguments is wrong!\n"
                << "[USAGE] ./regression_curve_5d <file_name>"
                << std::endl;
      return EXIT_FAILURE;
    }

    // ファイル名取得
    f_data = argv[1];

    // データ取得
    File file(f_data);
    if (!file.get_text(data)) {
      std::cout << "[ERROR] Failed to read the file!" << std::endl;
      return EXIT_FAILURE;
    }

    // データ一覧出力
    std::cout << std::fixed << std::setprecision(4);
    std::cout << "説明変数 X  目的変数 Y" << std::endl;
    for (i = 0; i < data.size(); i++)
      std::cout << std::setw(10) << std::right << data[i][0]
                << "  "
                << std::setw(10) << std::right << data[i][1]
                << std::endl;

    // 計算
    Calc calc(data);
    if (!calc.reg_curve_5d(a, b, c, d, e, f)) {
      std::cout << "[ERROR] Failed to calculate!" << std::endl;
      return EXIT_FAILURE;
    }

    // 結果出力
    std::cout << std::fixed << std::setprecision(8);
    std::cout << "---\n"
              << "a = " << std::setw(16) << std::right << a
              << "\n"
              << "b = " << std::setw(16) << std::right << b
              << "\n"
              << "c = " << std::setw(16) << std::right << c
              << "\n"
              << "d = " << std::setw(16) << std::right << d
              << "\n"
              << "e = " << std::setw(16) << std::right << e
              << "\n"
              << "f = " << std::setw(16) << std::right << f
              << std::endl;
  } catch (...) {
      std::cerr << "EXCEPTION!" << std::endl;
      return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}