mk-mode.com

mk-mode.com

Linux, Debian, IT, Server, PG, Ruby, Rails, Python, C++, Fortran, PC, MariaDB, math, GIS, etc...

C++ - ロジスティック回帰分析！

Updated:


少し前に、説明変数K個・目的変数1個のロジスティック回帰分析のアルゴリズムを Ruby で実装したことを紹介しました。

今回は、説明変数2個・目的変数1個のロジスティック回帰分析のアルゴリズムを C++ で実装してみました。

0. 前提条件

  • Debian GNU/Linux 11.5 (64bit) での作業を想定。
  • GCC 12.1.0 (G++ 12.1.0) (C++17) でのコンパイルを想定。

1. アルゴリズム

当ブログ過去記事を参照。

2. ソースコードの作成

  • ファイル読み込み部分、計算部分、実行部分とソースファイルを分けている。

File: file.hpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

#ifndef REGRESSION_LOGISTIC_FILE_HPP_
#define REGRESSION_LOGISTIC_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

File: file.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38

#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);
      }

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

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

File: calc.hpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

#ifndef REGRESSION_LOGISTIC_CALC_HPP_
#define REGRESSION_LOGISTIC_CALC_HPP_

#include <Eigen/Core>  // for 行列定義・行列演算
#include <vector>

class Calc {
  std::vector<std::vector<double>> data;                   // 元データ
  bool sigmoid(const Eigen::VectorXd&, Eigen::VectorXd&);  // シグモイド関数
  double cross_entropy_loss(
      const Eigen::VectorXd&, const Eigen::VectorXd&);     // 交差エントロピー誤差関数

public:
  Calc(std::vector<std::vector<double>>& data) : data(data) {}
                                            // コンストラクタ
  bool reg_logistic(std::vector<double>&);  // ロジスティック回帰（説明変数2個）の計算
};

#endif

File: calc.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134

#include "calc.hpp"

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

// 定数
static constexpr double       kAlpha = 0.01;      // 学習率
static constexpr double       kEps   = 1.0e-12;   // 閾値
static constexpr unsigned int kLoop  = 10000000;  // 最大ループ回数
static constexpr double       kBeta  = 5.0;       // 初期値: β
static constexpr double       kCel   = 0.0;       // 初期値: 交差エントロピー誤差
static constexpr double       kOne   = 1.0;       // 1

/**
 * @brief      ロジスティック回帰（説明変数2個）の計算
 *
 * @param[ref] パラメータ（定数・係数） ps (vector<double>)
 * @return     真偽(bool)
 * @retval     true  成功
 * @retval     false 失敗
 */
bool Calc::reg_logistic(std::vector<double>& ps) {
  std::size_t   e;  // 元の数
  std::size_t   n;  // サンプル数
  unsigned int  i;  // loop インデックス
  unsigned int  j;  // loop インデックス
  double     loss;  // 交差エントロピー誤差
  double loss_pre;  // 交差エントロピー誤差（退避用）

  try {
    // 元の数, サンプル数
    e = data[0].size() - 1;
    n = data.size();

    // データの行列化
    Eigen::MatrixXd mtx(n, 3); // n x 3 行列(double) b 用
    for (i = 0; i < n; i++)
      for (j = 0; j < 3; j++) mtx(i, j) = data[i][j];

    // β初期値 (e + 1 次元ベクトル)
    Eigen::VectorXd bs(e + 1);
    bs.setConstant(kBeta);

    // X の行列 (n 行 e + 1 列)
    // (第1列(x_0)は定数項なので 1 固定)
    Eigen::MatrixXd xs(n, e + 1);
    xs.setConstant(kOne);
    for (j = 0; j < e; j++) xs.col(j + 1) = mtx.col(j);

    // t のベクトル (n 次元ベクトル)
    Eigen::VectorXd ts(n);
    ts = mtx.col(e).transpose();

    //# 交差エントロピー誤差初期値
    loss = kCel;

    // y のベクトル (n 次元ベクトル)
    Eigen::VectorXd ys(n);

    // dE のベクトル (e + 1 次元ベクトル)
    Eigen::VectorXd des(e + 1);
    for (i = 0; i < kLoop; i++) {
      // シグモイド関数適用（予測値計算） (n 次元ベクトル)
      if (!sigmoid(xs * bs, ys)) {
        std::cout << "[ERROR] Failed to calculate! (in sigmoid() function)"
                  << std::endl;
        return EXIT_FAILURE;
      }
      // dE 計算 (e + 1 次元ベクトル)
      des = (ys - ts).transpose() * xs / n;
      // β 更新 (e + 1 次元ベクトル)
      bs -= kAlpha * des;
      // 前回算出交差エントロピー誤差退避
      loss_pre = loss;
      // 交差エントロピー誤差計算
      loss = cross_entropy_loss(ts, ys);
      // 今回と前回の交差エントロピー誤差の差が閾値以下になったら終了
      if (abs(loss - loss_pre)< kEps) break;
    }

    // 計算結果用 vector へ格納
    for (i = 0; i < e + 1; i++) ps.push_back(bs(i));
  } catch (...) {
    return false;  // 計算失敗
  }

  return true;  // 計算成功
}

/**
 * @brief      シグモイド関数
 *
 * @param[ref] 計算前ベクトル as (Eigen::VectorXd)
 * @param[ref] 計算後ベクトル ys (Eigen::VectorXd)
 * @return     真偽(bool)
 * @retval     true  成功
 * @retval     false 失敗
 */
bool Calc::sigmoid(const Eigen::VectorXd& as, Eigen::VectorXd& ys) {
  try {
    ys  = 1.0 / (1.0 + (-as).array().exp());
  } catch (...) {
    return false;  // 計算失敗
  }

  return true;  // 計算成功
}

/**
 * @brief      交差エントロピー誤差関数
 *
 * @param[ref] 実際の値(0 or 1)ベクトル        ts (Eigen::VectorXd)
 * @param[ref] 回帰式から得られる y 値ベクトル ys (Eigen::VectorXd)
 * @return     真偽(bool)
 * @retval     true  成功
 * @retval     false 失敗
 */
double Calc::cross_entropy_loss(
    const Eigen::VectorXd& ts, const Eigen::VectorXd& ys) {
  double loss;
  Eigen::VectorXd ones;  // all 1.0 のベクトル

  try {
    ones = Eigen::VectorXd::Ones(ys.rows());
    loss = (ts.array() * ys.array().log()
         + (ones - ts).array() * (ones - ys).array().log()).sum();
  } catch (...) {
    throw;
  }

  return loss;
}

File: regression_logistic.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73

/***********************************************************
  ロジスティック回帰計算
  * 説明(独立)変数2個、目的(従属)変数1個 限定

    DATE        AUTHOR       VERSION
    2022.11.02  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;  // データ配列
  unsigned int i;                         // loop インデックス
  std::vector<double> ps;                 // 定数・係数 beta

  try {
    // コマンドライン引数のチェック
    if (argc != 2) {
      std::cerr << "[ERROR] Number of arguments is wrong!\n"
                << "[USAGE] ./regression_logistic <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 << "説明変数 X  説明変数 Y  目的変数 Z" << std::endl;
    std::cout << std::fixed << std::setprecision(4);
    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_logistic(ps)) {
      std::cout << "[ERROR] Failed to calculate!" << std::endl;
      return EXIT_FAILURE;
    }

    // 結果出力
    std::cout << "betas =" << std::endl;
    std::cout << std::fixed << std::setprecision(8);
    for (auto p: ps)
      std::cout << std::setw(16) << std::right << p << std::endl;
  } catch (...) {
      std::cerr << "EXCEPTION!" << std::endl;
      return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

3. ソースコードのコンパイル

まず、以下のように Makefile を作成する。（行頭のインデントはタブ文字）

File: Makefile

gcc_options = -std=c++17 -Wall -O2 --pedantic-errors \
-I /usr/local/include/eigen-3.4.0

regression_logistic: regression_logistic.o file.o calc.o
  g++ $(gcc_options) -o $@ $^

regression_logistic.o : regression_logistic.cpp
  g++ $(gcc_options) -c $<

file.o : file.cpp
  g++ $(gcc_options) -c $<

calc.o : calc.cpp
  g++ $(gcc_options) -c $<

run : regression_logistic
  ./regression_logistic

clean :
  rm -f ./regression_logistic
  rm -f ./*.o

.PHONY : run clean

そして、ビルド（コンパイル＆リンク）。

$ make

4. 動作確認

まず、以下のような入力ファイルを用意する。
（各行は x, y （説明変数）と z （目的変数）の値）

File: data.txt

30 21 0
22 10 0
26 25 0
14 20 0
 6 10 1
 2 15 1
 6  5 1
10  5 1
19 15 1

そして、コマンドライン引数にデータファイルを指定して実行。

$ ./regression_logistic data.txt
説明変数 X  説明変数 Y  目的変数 Z
   30.0000     21.0000      0.0000
   22.0000     10.0000      0.0000
   26.0000     25.0000      0.0000
   14.0000     20.0000      0.0000
    6.0000     10.0000      1.0000
    2.0000     15.0000      1.0000
    6.0000      5.0000      1.0000
   10.0000      5.0000      1.0000
   19.0000     15.0000      1.0000
betas =
      8.99376352
     -0.30784038
     -0.26681572

betas の値が、上から \(\beta_0, \beta_1, \beta_2\) となっている。
また、前回の Ruby での実装時と同じ値になっていることも確認できた。

以上。





 

Sponsored Link

 

Comments

You May Also Enjoy

Debian 12 (bookworm) - Ruby 3.3 インストール（ソースビルド）！  (2024-05-18)

Debian GNU/Linux 12 (bookworm) に Ruby 3.3 をソースをビルドしてインストール方法についての記録です。 以前古いバージョンでの作業時に残していた記録を参考に作業を行い、今回更新した作業記録を貼付する形式の内容となっています。 （当然ながら、興味がなければスルーしてください）