Skip to content

Instantly share code, notes, and snippets.

@komasaru
Last active October 8, 2021 11:50
Show Gist options
  • Save komasaru/a10d55af6d91bddcda8a6d22a455d6ce to your computer and use it in GitHub Desktop.
Save komasaru/a10d55af6d91bddcda8a6d22a455d6ce to your computer and use it in GitHub Desktop.
C++ source code to compute multiple regression equations.(2d)
#include "calc.hpp"
#include <cmath>
#include <iostream>
#include <sstream>
#include <vector>
/**
* @brief 重回帰式(説明変数2個; 2次多項式モデル)の計算
*
* @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_multi_2e_2d(
double& a, double& b, double& c, double& d, double& e, double& f
) {
unsigned int i; // loop インデックス
double s_x1 = 0.0; // sum(x1 )
double s_x1x1 = 0.0; // sum(x1 * x1)
double s_x1x2 = 0.0; // sum(x1 * x2)
double s_x1x3 = 0.0; // sum(x1 * x3)
double s_x1x4 = 0.0; // sum(x1 * x4)
double s_x1x5 = 0.0; // sum(x1 * x5)
double s_x1y = 0.0; // sum(x1 * y )
double s_x2 = 0.0; // sum(x2 )
double s_x2x2 = 0.0; // sum(x2 * x2)
double s_x2x3 = 0.0; // sum(x2 * x3)
double s_x2x4 = 0.0; // sum(x2 * x4)
double s_x2x5 = 0.0; // sum(x2 * x5)
double s_x2y = 0.0; // sum(x2 * y )
double s_x3 = 0.0; // sum(x3 )
double s_x3x3 = 0.0; // sum(x3 * x3)
double s_x3x4 = 0.0; // sum(x3 * x4)
double s_x3x5 = 0.0; // sum(x3 * x5)
double s_x3y = 0.0; // sum(x3 * y )
double s_x4 = 0.0; // sum(x4 )
double s_x4x4 = 0.0; // sum(x4 * x4)
double s_x4x5 = 0.0; // sum(x4 * x5)
double s_x4y = 0.0; // sum(x4 * y )
double s_x5 = 0.0; // sum(x5 )
double s_x5x5 = 0.0; // sum(x5 * x5)
double s_x5y = 0.0; // sum(x5 * y )
double s_y = 0.0; // sum(y )
double x1 = 0.0; // x1 計算用
double x2 = 0.0; // x2 計算用
double x3 = 0.0; // x3 計算用
double x4 = 0.0; // x4 計算用
double x5 = 0.0; // x5 計算用
double y = 0.0; // y 計算用
try {
// データ数
cnt = data.size();
// sum(x1), sum(x1 * x1), sum(x1 * x2), ...
for (i = 0; i < cnt; i++) {
x1 = data[i][0];
x2 = data[i][1];
x3 = x1 * x2;
x4 = x1 * x1;
x5 = x2 * x2;
y = data[i][2];
s_x1 += x1;
s_x1x1 += x1 * x1;
s_x1x2 += x1 * x2;
s_x1x3 += x1 * x3;
s_x1x4 += x1 * x4;
s_x1x5 += x1 * x5;
s_x1y += x1 * y;
s_x2 += x2;
s_x2x2 += x2 * x2;
s_x2x3 += x2 * x3;
s_x2x4 += x2 * x4;
s_x2x5 += x2 * x5;
s_x2y += x2 * y;
s_x3 += x3;
s_x3x3 += x3 * x3;
s_x3x4 += x3 * x4;
s_x3x5 += x3 * x5;
s_x3y += x3 * y;
s_x4 += x4;
s_x4x4 += x4 * x4;
s_x4x5 += x4 * x5;
s_x4y += x4 * y;
s_x5 += x5;
s_x5x5 += x5 * x5;
s_x5y += x5 * y;
s_y += y;
}
// 行列1行目
mtx.push_back({(double)cnt, s_x1, s_x2, s_x3, s_x4, s_x5, s_y});
// 行列2行目
mtx.push_back({mtx[0][1], s_x1x1, s_x1x2, s_x1x3, s_x1x4, s_x1x5, s_x1y});
// 行列3行目
mtx.push_back({mtx[0][2], mtx[1][2], s_x2x2, s_x2x3, s_x2x4, s_x2x5, s_x2y});
// 行列4行目
mtx.push_back({mtx[0][3], mtx[1][3], mtx[2][3], s_x3x3, s_x3x4, s_x3x5, s_x3y});
// 行列5行目
mtx.push_back({mtx[0][4], mtx[1][4], mtx[2][4], mtx[3][4], s_x4x4, s_x4x5, s_x4y});
// 行列6行目
mtx.push_back({mtx[0][5], mtx[1][5], mtx[2][5], mtx[3][5], mtx[4][5], s_x5x5, s_x5y});
// 計算(ガウスの消去法)
if (!solve_ge(mtx)) {
std::cout << "[ERROR] Failed to solve by the Gauss-Ellimination method!"
<< std::endl;
return false;
}
// b0, ..., b5
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 (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++) {
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; // 計算成功
}
#ifndef REGRESSION_MULTI_2E_2D_CALC_HPP_
#define REGRESSION_MULTI_2E_2D_CALC_HPP_
#include <vector>
class Calc {
std::vector<std::vector<double>> data; // 元データ
std::vector<std::vector<double>> mtx; // 計算用行列
bool solve_ge(std::vector<std::vector<double>>&); // ガウスの消去法
public:
Calc(std::vector<std::vector<double>>& data) : data(data) {}
unsigned int cnt; // データ件数
bool reg_multi_2e_2d(double&, double&, double&, double&, double&, double&);
// 重回帰式(説明変数2個; 2次多項式モデル)の計算
};
#endif
#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, z;
while (iss >> x >> y >> z) {
rec.push_back(x);
rec.push_back(y);
rec.push_back(z);
}
// 1行分ベクタを data ベクタに追加
if (rec.size() != 0) data.push_back(rec);
}
return true;
} catch (...) {
std::cerr << "EXCEPTION!" << std::endl;
return false;
}
return true; // 読み込み成功
}
#ifndef REGRESSION_MULTI_2E_2D_FILE_HPP_
#define REGRESSION_MULTI_2E_2D_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
/***********************************************************
重回帰式計算(説明(独立)変数2個、2次多項式モデル)
* y = b0 + b1x1 + b2x2 + b3x1x2 + b4x1^2 + b5x2^2
* y = b0 + b1x1 + b2x2 + b3x3 + b4x4 + b5x5
(x3 = x1x2, x4 = x1^2, x5 = x2^2)
ということ。
DATE AUTHOR VERSION
2020.07.15 mk-mode.com 1.00 新規作成
Copyright(C) 2020 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; // 定数 b0
double b; // 係数 b1
double c; // 係数 b2
double d; // 係数 b3
double e; // 係数 b4
double f; // 係数 b5
try {
// コマンドライン引数のチェック
if (argc != 2) {
std::cerr << "[ERROR] Number of arguments is wrong!\n"
<< "[USAGE] ./regression_multi_2e_2d <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 目的変数 Z" << 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::setw(10) << std::right << data[i][2]
<< std::endl;
// 計算
Calc calc(data);
if (!calc.reg_multi_2e_2d(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"
<< "b0 = " << std::setw(16) << std::right << a
<< "\n"
<< "b1 = " << std::setw(16) << std::right << b
<< "\n"
<< "b2 = " << std::setw(16) << std::right << c
<< "\n"
<< "b3 = " << std::setw(16) << std::right << d
<< "\n"
<< "b4 = " << std::setw(16) << std::right << e
<< "\n"
<< "b5 = " << std::setw(16) << std::right << f
<< std::endl;
} catch (...) {
std::cerr << "EXCEPTION!" << std::endl;
return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment