C言語による重回帰分析の実装：最小二乗法を用いた多変量データ解析について解説

この記事ではC言語を用い、最小二乗法により重回帰分析を実装する手順を解説します。

多変量データから、データ行列$X$と目的変数$y$をもとに推定係数$\hat{\beta }=(X^{T}X{)}^{-1}{X}^{T}y$を計算する流れを紹介し、プログラムの作成方法について具体例を交えて説明します。

重回帰分析の基本

多変量データの特徴

多変量データとは、複数の変数を同時に扱うデータであり、各観測対象ごとに複数の値が対応しているデータ形式です。

例えば、住宅価格の予測では、家の広さ、築年数、立地など複数の要因が影響するため、各要因が列として管理されます。

そのため、各変数間の相互関係を考慮し、全体としての関係性をモデル化することが求められます。

最小二乗法の原理

最小二乗法は、予測値と実際の値の誤差の二乗和を最小化することを目的として、回帰係数を求める手法です。

具体的には、実測値 $y_i$ と予測値 $\widehat{y_i}$ の差の二乗和

$$\epsilon =\sum _{i=1}^n(y_i–\widehat{y_i}{)}^2$$

を最小にする回帰係数 $\mathit{\beta }$ を求めます。

この手法は、ノイズが含まれる実データに対しても、モデルのフィッティングを比較的安定して行える点が魅力です。

正規方程式の導出

最小二乗誤差を最小化するために、偏微分を用いて誤差関数の最小値を求めると、正規方程式と呼ばれる

$${\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\mathit{\beta }={\mathbf{X}}^{\mathrm{\top }}\mathbf{y}$$

が得られます。

ここで、$\mathbf{X}$ はデザイン行列、$\mathbf{y}$ は目的変数のベクトルです。

この方程式を解くことで、重回帰分析に必要な係数 $\mathit{\beta }$ を計算することができます。

C言語による実装概要

プログラム構成と全体の流れ

C言語での重回帰分析の実装は、以下のような流れで構成されます。

1. データの入力および前処理を行い、デザイン行列 $\mathbf{X}$ と目的変数ベクトル $\mathbf{y}$ を構築する。
2. 行列演算のための各種関数(転置、積、逆行列計算)を利用して、正規方程式 ${\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\mathit{\beta }={\mathbf{X}}^{\mathrm{\top }}\mathbf{y}$ を解く。
3. 計算結果として得られた回帰係数 $\mathit{\beta }$ を出力する。

全体の流れとしては、モジュールごとに機能を分割し、行列計算部分とデータ入出力部分を明確に分けることで、コードの可読性と保守性を向上させています。

使用するライブラリと主要な関数

この実装では、標準ライブラリを主に利用します。

対象となるライブラリと関数は以下の通りです。

• ライブラリ
  • stdio.h : 入出力処理のために使用
  • stdlib.h : 動的メモリ確保や標準ユーティリティ関数のために使用
  • math.h : 数学的計算(場合によっては平方根やその他の関数)のために使用
• 主要な関数
  • readData() : データの入力および行列の構築を行う関数
  • matrixTranspose() : 行列の転置を計算する関数
  • matrixMultiply() : 行列の積を計算する関数
  • matrixInverse() : 逆行列を計算する関数(ガウス・ジョルダン法などを内部で利用)
  • computeBeta() : 正規方程式を構築し、回帰係数 $\mathit{\beta }$ を求める関数
  • printResults() : 計算結果の出力を行う関数

行列演算の実装詳細

行列の転置と積の計算

行列の転置処理

行列の転置処理では、元の行列の各要素を行と列を入れ替えた位置にコピーします。

以下は、シンプルな行列転置のサンプルコードです。

#include <stdio.h>
#define ROWS 2
#define COLS 3
// matrixTranspose : 行列の転置を実施する関数
void matrixTranspose(double input[ROWS][COLS], double output[COLS][ROWS]) {
    int i, j;
    for (i = 0; i < ROWS; i++) {
        for (j = 0; j < COLS; j++) {
            output[j][i] = input[i][j];
        }
    }
}
int main(void) {
    double A[ROWS][COLS] = {
        {1.0, 2.0, 3.0},
        {4.0, 5.0, 6.0}
    };
    double AT[COLS][ROWS]; // 転置後の行列
    matrixTranspose(A, AT);
    // 転置行列の表示
    int i, j;
    printf("Transpose of matrix A:\n");
    for (i = 0; i < COLS; i++) {
        for (j = 0; j < ROWS; j++) {
            printf("%lf ", AT[i][j]);
        }
        printf("\n");
    }
    return 0;
}

Transpose of matrix A:
1.000000 4.000000
2.000000 5.000000
3.000000 6.000000

行列積のアルゴリズム

行列積は、第一行列の各行と第二行列の各列の内積を求めることで計算します。

アルゴリズムは、外側の二重ループで結果行列の行と列を走査し、内側のループで対応する要素を掛け合わせて足し合わせます。

以下のサンプルコードは、二つの行列の積を計算する方法を示しています。

#include <stdio.h>
#define M 2   // Aの行数
#define N 3   // Aの列数、Bの行数
#define P 2   // Bの列数
// matrixMultiply : 行列AとBの積を計算する関数
void matrixMultiply(double A[M][N], double B[N][P], double C[M][P]) {
    int i, j, k;
    for (i = 0; i < M; i++) {
        for (j = 0; j < P; j++) {
            C[i][j] = 0.0;
            for (k = 0; k < N; k++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }
}
int main(void) {
    double A[M][N] = {
        {1.0, 2.0, 3.0},
        {4.0, 5.0, 6.0}
    };
    double B[N][P] = {
        {7.0, 8.0},
        {9.0, 10.0},
        {11.0, 12.0}
    };
    double C[M][P];  // 積の結果の行列
    matrixMultiply(A, B, C);
    // 結果行列の表示
    int i, j;
    printf("Result of A * B:\n");
    for (i = 0; i < M; i++) {
        for (j = 0; j < P; j++) {
            printf("%lf ", C[i][j]);
        }
        printf("\n");
    }
    return 0;
}

Result of A * B:
58.000000 64.000000
139.000000 154.000000

逆行列の計算

ガウス・ジョルダン法の概要

逆行列の計算では、ガウス・ジョルダン法を用いて行列を単位行列に変換する操作を行います。

具体的には、元の行列に単位行列を連結し、前進消去と後退代入を組み合わせることで、元の行列を単位行列に変換すると同時に、連結している部分が逆行列に変換されます。

この方法は、行列が特異(非可逆)な場合には適用できないため、計算前に行列の性質を確認する必要があります。

計算上の留意点

ガウス・ジョルダン法を実装する際は、以下の点に注意が必要です。

• ピボット要素がゼロまたは極小の場合、数値的安定性が損なわれるため、行の入れ替えなどを行う工夫が必要です。
• 浮動小数点演算の誤差により、厳密なゼロ判定が難しいため、十分小さい値との比較を行います。
• 行列が特異である場合、逆行列の計算はできないため、エラー処理を適切に実装します。

重回帰分析計算手順

データ入力と行列の構築

重回帰分析では、まずデータを入力し、デザイン行列 $\mathbf{X}$ と目的変数ベクトル $\mathbf{y}$ を構築します。

各サンプルごとに、説明変数を行ベクトルとしてまとめ、バイアス項を追加することで、$\mathbf{X}$ は

$$\mathbf{X}=\left[\begin{array}{cccc}1& x_{11}& x_{12}& \cdots \\ 1& x_{21}& x_{22}& \cdots \\ ⋮& ⋮& ⋮& \ddots \end{array}\right]$$

という形になります。

データの入力は標準入力やファイルから行うことができ、要素の整合性が保たれるように注意が必要です。

推定係数計算の実装

正規方程式の構築と解法

推定係数 $\mathit{\beta }$ は、正規方程式

$${\mathbf{X}}^{\mathrm{\top }}\mathbf{X}\mathit{\beta }={\mathbf{X}}^{\mathrm{\top }}\mathbf{y}$$

を解くことで求めます。

計算手順は以下の通りです。

1. 行列 $\mathbf{X}$ の転置 ${\mathbf{X}}^{\mathrm{\top }}$ を求める。
2. 計算した ${\mathbf{X}}^{\mathrm{\top }}$ を用い、${\mathbf{X}}^{\mathrm{\top }}\mathbf{X}$ を算出する。
3. ${\mathbf{X}}^{\mathrm{\top }}\mathbf{X}$ の逆行列を求める。
4. 最終的に、$\mathit{\beta }=({\mathbf{X}}^{\mathrm{\top }}\mathbf{X}{)}^{-1}{\mathbf{X}}^{\mathrm{\top }}\mathbf{y}$ を計算する。

ここでは、既に実装している matrixTranspose()、matrixMultiply()、matrixInverse() などの関数を組み合わせて解法を実現します。

結果の出力処理

計算された回帰係数 $\mathit{\beta }$ を、整形して出力します。

出力処理では、有効数字の調整やフォーマットの整形を行い、ユーザーが結果を容易に確認できるようにします。

以下のサンプルコードは、回帰係数を出力する際の基本的な例です。

#include <stdio.h>
#include <stdlib.h>
// printResults : 回帰係数を表示する関数
void printResults(double beta[], int size) {
    int i;
    printf("Estimated coefficients (beta):\n");
    for (i = 0; i < size; i++) {
        printf("beta[%d] = %lf\n", i, beta[i]);
    }
}
int main(void) {
    // 仮の回帰係数データ
    double beta[] = {1.234, 2.345, 3.456};
    int size = sizeof(beta) / sizeof(beta[0]);
    printResults(beta, size);
    return 0;
}

Estimated coefficients (beta):
beta[0] = 1.234000
beta[1] = 2.345000
beta[2] = 3.456000

プログラム実行とデバッグ

実行例の紹介

プログラムの実行例として、予め用意したサンプルデータで重回帰分析を行った場合の結果を紹介します。

実行時には、データの読み込み、行列演算、そして回帰係数の計算が順次行われ、最終的に各係数が標準出力に表示されます。

例えば、以下のような出力が得られる場合があります。

• 入力例:

・説明変数とバイアス項を持つデザイン行列

・目的変数のベクトル

• 出力例:

・beta[0] = 0.5

・beta[1] = 2.3

・beta[2] = -1.7

これにより、各変数が目的変数にどのような影響を与えているかが定量的に理解できます。

エラーチェックとデバッグ方法

プログラム実行時に考慮すべきエラーとしては、以下の点があります。

• 行列のサイズや次元の不一致

入力データが正しく整形されていない場合は、計算が正しく行われず、エラーメッセージを出力する仕組みが必要です。

• 逆行列計算時の特異行列

${\mathbf{X}}^{\mathrm{\top }}\mathbf{X}$ が特異な場合、逆行列は存在しないため、エラー判定を行う必要があります。

この場合、適切なエラーメッセージを出力し、処理を中断するなどの対策が求められます。

• 数値誤差の管理

浮動小数点計算に伴う誤差を考慮し、十分小さい値との比較により、実質的なゼロ判定を行う手法を採用することが重要です。

デバッグの際には、各関数ごとにテスト用のデータを与えて正しいアルゴリズムが実装されているかを確認し、標準出力へのデバッグ用のログを活用するなど、段階的な検証が効果的です。

まとめ

この記事では、C言語を用いた重回帰分析の基礎概念、つまり多変量データの扱いや最小二乗法の原理、正規方程式の導出方法について説明しました。

また、行列の転置や積、ガウス・ジョルダン法による逆行列計算といった行列演算の実装詳細、デザイン行列の構築、回帰係数の算出手法、データ入出力の流れ、実行例およびエラー対策についても解説し、実践的なプログラム作成に役立つ内容となっています。