最小二乗法 #C#

最小二乗法 (least squares method) とは、測定で得られた数値の組を、
適当なモデルから想定される1次関数、対数曲線など特定の関数を用いて近似するときに、
想定する関数が測定値に対してよい近似となるように、
残差の二乗和を最小とするような係数を決定する方法、あるいはそのような方法によって近似を行うことである
(引用 Wikipedia)

線形最小二乗法

XY 座標を表す構造体を定義

struct Coordinate
{
    public double X { get; set; }
    public double Y { get; set; }
}

一次方程式 (y = Ax + B) を表す構造体を定義

// y = Ax + B
struct LinearEquation
{
    public double A { get; set; }
    public double B { get; set; }

    public double Solve(double x) => this.A * x + this.B;
}

LinearEquation クラスに最小二乗法の静的メソッドを追加

public static LinearEquation Approximate(Coordinate[] data)
{
    if (data == null) throw new ArgumentNullException(nameof(data));
    var n = data.Length;
    double Σx = 0, Σy = 0, Σxx = 0, Σxy = 0;
    for (int i = 0; i < n; i++)
    {
        var x = data[i].X;
        var y = data[i].Y;
        Σx += x;
        Σy += y;
        Σxx += Math.Pow(x, 2);
        Σxy += x * y;
    }
    return new LinearEquation()
    {
        A = (n * Σxy - Σx * Σy) / (n * Σxx - Math.Pow(Σx, 2)),
        B = (Σy - Σx * a) / n,
    };
}

最小二乗法の手順

フィッティングさせたい数式を決定 (例では y = Ax + B)
二乗和の式を展開
各係数で偏微分
連立方程式を解いて係数の値を決定

導出例

線形方程式を定義

f_x = Ax + B \\
\tag{1}

二乗和の式に式 (1) を代入して展開

\begin{align}
\sum_{i=0}^{n-1} (y_i-f_i)^2
&= \sum{}y^2 - \sum{}2yf_x + \sum{}f_x^2 \\
&= \sum{}y^2 - 2A\sum{}xy - 2B\sum{}y + A^2\sum{}x^2 + 2AB\sum{}x + nB^2 \\
\end{align}
\tag{2}

式 (2) について、A と B のそれぞれで偏微分すると

\begin{align}
-2\sum{}xy + 2A\sum{}x^2 + 2B\sum{}x = 0 \\
-2\sum{}y + 2A\sum{}x + 2nB = 0 \\
\end{align}
\tag{3}

式 (3) を連立方程式として解くと

\begin{align}
A &= \frac{n\sum{}xy - \sum{}x \sum{}y}{n\sum{}x^2 - (\sum{}x)^2} \\
B &= \frac{\sum{}y - a\sum{}x}{n} \\
\end{align}
\tag{4}

指数関数に応用

指数方程式 (2-1) について

y = B･e^{Ax}
\tag{2-1}

両辺の対数をとって (2-2) の形に変形

log_e{y} = log_e{B} + Ax \\
\tag{2-2}

対数部分を式 (2-3) の変数に置き換え

\begin{align}
y\,' &= log_e{y} \\
B\,' &= log_e{B} \\
\end{align}
\tag{2-3}

一次方程式と同じ式 (2-4) となる

y\,' = Ax + B\,' \\
\tag{2-4}

これをコードで記述すると、下記のようにシンプルに！

// y = B･exp(Ax)
struct ExponentialEquation
{
    public double A { get; set; }
    public double B { get; set; }

    public double Solve(double x) => Math.Exp(this.A * x) * this.B;

    public static ExponentialEquation Approximate(Coordinate[] data)
    {
        var pt = data.Select(item => new Coordinate(item.X, Math.Log(item.Y)));
        var eq = LinearEquation.Approximate(pt.ToArray());
        return new ExponentialEquation() { A = eq.A, B = Math.Exp(eq.B) };
    }
}