C++で与えられた点集合を最良の直線で表現する方法（最小二乗法）

本記事では、点の集合を最も適切な直線で表現する方法について解説します。点の集合の座標 (x, y) が与えられたとき、最良の直線 y = mx + c を求めることがゴールです。つまり、傾き m と切片 c の値を見つけるだけです。まずは例を見てみましょう。

入力: no_of_points = 4
x1 = 2, y1 = 3,
x2 = 5, y2 = 6,
x3 = 1, y3 = 3,
x4 = 4, y4 = 5

出力: m = 0.8, c = 1.85
説明: 求めた m と c の値を任意の点 (xi, yi) に対する方程式 y = mx + c に代入すると、すべての点をできるだけよく通る直線が得られます。
(x2, y2) に m と c を代入すると、
左辺: mx + c = 0.8 × 5 + 1.85 = 5.85
右辺: y = 6 となり、左辺とほぼ一致します。

入力: no_of_points = 3
x1 = 3, y1 = 6,
x2 = 2, y2 = 4,
x3 = 1, y3 = 3

出力: m = 1.5, c = 1.33

解法のアプローチ

この問題を解くには、m と c の値を求める必要があります。点の数が 2 つの場合は一意な解が必ず存在しますが、点の数が 3 つ以上になると、すべての点を通る直線が存在するとは限りません。

点の数を n とすると、n 個の方程式 fn = mxn + c が得られます。この方程式が最良の近似となるためには、fi の値がそれぞれの yi と等しい、またはできるだけ近い値になるようにする必要があります。

そこで Z = (fi − yi)² とおきます。すべての点についてこの値を最小化することを目指します。(fi − yi) を2乗するのは、正と負の誤差が打ち消し合わないようにするためです。

Z を最小にするためには、次の条件を満たす必要があります。

∂Z/∂m = 0 および ∂Z/∂c = 0

これらの方程式を解くと、次の 2 つの式が得られます。

Σy = m × Σx + n × c
Σxy = m × Σx² + c × Σx

これを m と c について整理すると、

m = (n × Σxy − Σx × Σy) / (n × Σx² − (Σx)²)
c = (Σy − m × Σx) / n

これで、最終的な方程式の m と c を直接求める公式が得られました。この手法は統計学では「最小二乗法」と呼ばれ、線形回帰の基礎となる重要な考え方です。

実装例

上記のアプローチをC++で実装したコード

#include <cmath>
#include <iostream>
using namespace std;
int main(){
    int X[] = { 3, 2, 1 };
    int Y[] = { 6, 4, 3};
    int no_of_points = sizeof(X) / sizeof(X[0]);
    float m, c;
    int sum_of_X = 0, sum_of_X2 = 0, sum_of_Y = 0, sum_of_XY = 0;
    // 方程式の各項を計算する
    for (int i = 0; i < no_of_points; i++) {
        sum_of_X = sum_of_X + X[i];
        sum_of_X2 = sum_of_X2 + pow(X[i],2);
        sum_of_Y = sum_of_Y + Y[i];
        sum_of_XY = sum_of_XY + (X[i] * Y[i]);
    }
    // 公式を使って m と c の値を計算する
    m = (no_of_points * sum_of_XY - sum_of_X * sum_of_Y) / (no_of_points * sum_of_X2 - pow(sum_of_X,2));
    c = (sum_of_Y - m * sum_of_X) / no_of_points;
    cout << "m = " << m;
    cout << "\nc = " << c;
    return 0;
}

出力

m = 1.5
c = 1.33333

まとめ

このチュートリアルでは、与えられた点の集合を表す最良近似直線を求める方法について解説しました。まず m と c の公式を導出し、それをそのまま適用するというシンプルなアプローチを紹介しました。また、この問題を解く C++ プログラムも示しましたが、同様の処理は C、Java、Python などの他のプログラミング言語でも簡単に実装できます。このチュートリアルがお役に立てば幸いです。