最小二乗法

データ群の関係性の関係式を導く方法だよ。

定義だよ

線形近似（<math>y=ax+b</math>）の場合だよ

<math>a=\displaystyle{\frac{n\sum_{i=1}^nx_iy_i-\sum_{i=1}^nx_i\sum_{i=1}^ny_i}{n\sum_{i=1}^n{x_i}^2-(\sum_{i=1}^nx_i)^2}}</math>

<math>b=\displaystyle{\frac{\sum_{i=1}^n{x_i}^2\sum_{i=1}^ny_i-\sum_{i=1}^nx_iy_i\sum_{i=1}^nx_i}{n\sum_{i=1}^n{x_i}^2-(\sum_{i=1}^nx_i)^2}}</math>

実は…

この計算式って、以下の行列を展開して計算した式だよ。
<math>\begin{pmatrix}\sum {x_i}^2 & \sum x_i \\\sum x_i & n \end{pmatrix} \begin{pmatrix}a \\b \end{pmatrix} = \begin{pmatrix}\sum x_iy_i\\\sum y_i \end{pmatrix}</math>
この計算式を解くには「ガウスの消去法」が分かりやすいんだよ。

特性だよ

残差の二乗和を最少にしているよ

残差って近似式からのyの値と実際のyの値の事の差分のことなんだ。それぞれのデータについて残差を求めて、その合計が最小になるように引いた直線が線形近似直線ってことになるよ。
二次曲線で最小になるように引いた場合は二次近似曲線だよ。

プログラミングだよ

エクセルだよ

線形近似だよ

傾きを求める場合だよ

• SLOPE([Ｙ軸の範囲],[Ｘ軸の範囲]) : 線形近似の傾きを求めるよ

切片を求める場合だよ

• INTERCEPT([Ｙ軸の範囲],[Ｘ軸の範囲]) : 線形近似の切片を求めるよ

多項近似だよ

例えば<math>y=ax^2+bx+c</math>の場合だよ
1.準備

y の範囲 A2:A10

x の範囲 B2:B10

とした場合…

c列にb列の二乗を算出（C2に「=B2^2」を入力しC10までコピー）するよ

2.算出

係数 a =INDEX(LINEST(A2:A10,B2:C10),1)

係数 b =INDEX(LINEST(A2:A10,B2:C10),2)

定数 c =INDEX(LINEST(A2:A10,B2:C10),3)

Ｃ＃だよ

線形近似だよ

計算式から算出に必要な要素は以下のモノだよ。
<math>\displaystyle{\sum_{i=1}^nx_iy_i}</math>、 <math>\displaystyle{\sum_{i=1}^nx_i}</math>、 <math>\displaystyle{\sum_{i=1}^ny_i}</math>、 <math>\displaystyle{\sum_{i=1}^n{x_i}^2}</math>
それぞれを算出した後に、傾きおよび切片を算出するよ。

//傾きを求める関数
private Double SLOPE(Int32[] p_IntArrayX, Int32[] p_IntArrayY)
{
    if (p_IntArrayX.Length != p_IntArrayY.Length)
    {
        throw new Exception();
    }

    Double XYSum = 0;
    Double XSum = 0;
    Double YSum = 0;
    Double XPowSum = 0;

    for (Int32 i = 0; i < p_IntArrayX.Length; i++)
    {
        XYSum += p_IntArrayX[i] * p_IntArrayY[i];
        XSum += p_IntArrayX[i];
        YSum += p_IntArrayY[i];
        XPowSum += p_IntArrayX[i] * p_IntArrayX[i];
    }

    return (p_IntArrayX.Length * XYSum - XSum * YSum) / (p_IntArrayX.Length * XPowSum - XSum * XSum);
}

//切片を求める関数
private Double INTERCEPT(Int32[] p_IntArrayX, Int32[] p_IntArrayY)
{
    if (p_IntArrayX.Length != p_IntArrayY.Length)
    {
        throw new Exception();
    }

    Double XYSum = 0;
    Double XSum = 0;
    Double YSum = 0;
    Double XPowSum = 0;

    for (Int32 i = 0; i < p_IntArrayX.Length; i++)
    {
        XYSum += p_IntArrayX[i] * p_IntArrayY[i];
        XSum += p_IntArrayX[i];
        YSum += p_IntArrayY[i];
        XPowSum += p_IntArrayX[i] * p_IntArrayX[i];
    }

    return (XPowSum * YSum - XYSum * XSum)/(p_IntArrayX.Length * XPowSum - XSum * XSum);
}

多項近似だよ

毎回毎回計算式を導出してプログラムを組むのは大変だよね…
だから、多項近似ができるプログラムを上げておくよ。
使うときに注意点があって…変数の範囲に気を付けてね
必要に応じて変数の型を変更すると良いと思うよ

// 係数算出処理だよ
public Double[] CoefficientCalc(List<InputElement> p_InputElementList,Int32 p_Dimension)
{
    //求める係数の数は次数＋１だよね
    p_Dimension += 1;

    // ガウスの消去法で解く配列の作成をするよ
    Double[,] l_A = new Double[p_Dimension, p_Dimension + 1];
    for (Int32 i = 0; i < p_Dimension; i++)
    {
        for (Int32 j = 0; j < p_Dimension + 1; j++)
        {
            l_A[i,j] = new Double();
        }
    }

    // 値の格納をするよ
    for (Int32 i = 0; i < p_Dimension; i++)
    {
        for (Int32 j = 0; j < p_Dimension; j++)
        {
            for (Int32 k = 0; k < p_InputElementList.Count; k++)
            {
                l_A[i, j] += Math.Pow(p_InputElementList[k].XElement, i + j);
            }
        }
    }
    for (Int32 i = 0; i < p_Dimension; i++)
    {
        for (Int32 k = 0; k < p_InputElementList.Count; k++)
        {
            l_A[i, p_Dimension] += Math.Pow(p_InputElementList[k].XElement, i) * p_InputElementList[k].YElement;
        }
    }


    return this.Gauss(l_A, p_Dimension);
}

// ガウスの消去法処理だよ
private Double[] Gauss(Double[,] p_A, Int32 p_Dimension)
{
    //係数でソートするよ
    for (Int32 i = 0; i < p_Dimension; i++)
    {
        Double l_m = 0;
        Int32 l_pivot = i;

        //係数が一番大きい行を探すよ
        for (Int32 l = i; l < p_Dimension; l++)
        {
            if (l_m < Math.Abs(p_A[l, i]))
            {
                l_m = Math.Abs(p_A[l, i]);
                l_pivot = l;
            }
        }

        //行の入れ替えするよ
        if (l_pivot != i)
        {
            Double l_b = 0;
            for (Int32 j = 0; j < p_Dimension + 1; j++)
            {
                l_b = p_A[i, j];
                p_A[i, j] = p_A[l_pivot, j];
                p_A[l_pivot, j] = l_b;
            }
        }
    }

    // 前進消去処理だよ
    for (Int32 k = 0; k < p_Dimension; k++)
    {
        Double l_p = p_A[k, k];
        p_A[k, k] = 1;

        for (Int32 j = k + 1; j < p_Dimension + 1; j++)
        {
            p_A[k, j] /= l_p;
        }

        for (Int32 i = k + 1; i < p_Dimension; i++)
        {
            Double l_q = p_A[i, k];

            for (Int32 j = k + 1; j < p_Dimension + 1; j++)
            {
                p_A[i, j] -= l_q * p_A[k, j];
            }

            p_A[i, k] = 0;
        }
    }

    // 後退代入処理だよ
    Double[] l_Result = new Double[p_Dimension];
    for (Int32 i = 0; i < p_Dimension; i++)
    {
        l_Result[i] = new Double();
    }

    for (Int32 i = p_Dimension - 1; 0 <= i; i--)
    {
        l_Result[i] = p_A[i, p_Dimension];
        for (Int32 j = p_Dimension - 1; i < j; j--)
        {
            l_Result[i] -= p_A[i, j] * l_Result[j];
        }
    }

    return l_Result;
}

// データ投入用のクラスだよ
public class InputElement
{
    private Double m_XElement;
    private Double m_YElement;

    public InputElement()
    {
        this.m_XElement = 0;
        this.m_YElement = 0;
    }

    public InputElement(Double p_XElement, Double p_YElement)
    {
        this.m_XElement = p_XElement;
        this.m_YElement = p_YElement;
    }

    public Double XElement
    {
        get
        {
            return this.m_XElement;
        }
        set
        {
            this.m_XElement = value;
        }
    }

    public Double YElement
    {
        get
        {
            return this.m_YElement;
        }
        set
        {
            this.m_YElement = value;
        }
    }
}

参考サイト

最小二乗法の数値計算例（n次多項式で近似）（Ｃ言語）