直線フィッティング - わかばテクノロジ

最適化
直線フィッティング; 線形回帰; Linear Regression

概要

直線フィッティングとは独立変数\(x_i\)と従属変数\(y_i\)の組で表される\(N\)個のデータ\(i = 1 \cdots N \)があるとき、
直線\(\tilde{y_i} = a + b x_i\)または原点を通る直線\(\tilde{y_i} = a x_i\)と、\(y_i\)との誤差の総和を最小化する、\(a, b\)の組または\(a\)を求めることである。
誤差には2乗和が用いられる。
\(\quad \displaystyle E = \sum_{i=1}^{N} (y_i - \tilde{y_i})^2 \)

誤差最小値における\(a, b\)の導出

最小化する対象である誤差の2乗和は以下である。
\(\begin{eqnarray} \quad \displaystyle E &=& \sum_{i=1}^{N} (y_i - a - b x_i)^2 \\ \quad &=& \sum_{i=1}^{N} y_i^2-2 b x_i y_i-2 a y_i+b^2 x_i^2+2 a b x_i+a^2 \\ \quad &=& \sum_{i=1}^{N} y_i^2-2 b \sum_{i=1}^{N} x_i y_i-2 a \sum_{i=1}^{N} y_i+b^2 \sum_{i=1}^{N} x_i^2+2 a b \sum_{i=1}^{N} x_i+ a^2 N \\ \end{eqnarray} \)
\(E\)が最小化であるとき、\(a, b\)について偏微分はともに0になる。
\(\begin{eqnarray} \quad \displaystyle \frac{\partial E}{\partial a} &=& -2 \sum_{i=1}^{N} y_i+2 b \sum_{i=1}^{N} x_i+2 a N = 0 \quad &\rightarrow& a N + b \sum_{i=1}^{N} x_i = \sum_{i=1}^{N} y_i \\ \quad \displaystyle \frac{\partial E}{\partial b} &=& -2 \sum_{i=1}^{N} x_i y_i + 2 b \sum_{i=1}^{N} x_i^2 + 2 a \sum_{i=1}^{N} x_i = 0 \quad &\rightarrow& a \sum_{i=1}^{N} x_i + b \sum_{i=1}^{N} x_i^2 = \sum_{i=1}^{N} x_i y_i \\ \end{eqnarray}\)
2式を連立させて\(a, b\)について解く。
\(\quad \begin{pmatrix} N & \sum_{i=1}^{N} x_i \\ \sum_{i=1}^{N} x_i & \sum_{i=1}^{N} x_i^2 \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} \sum_{i=1}^{N} y_i \\ \sum_{i=1}^{N} x_i y_i \end{pmatrix} \)
であるから
\( \begin{eqnarray} \quad a &=& \frac{ \displaystyle \sum_{i=1}^{N} x_i \sum_{i=1}^{N} x_i y_i - \sum_{i=1}^{N} x_i^2 \sum_{i=1}^{N} y_i }{ \left( \displaystyle \sum_{i=1}^{N} x_i \right)^2 - N \displaystyle \sum_{i=1}^{N} x_i^2 } \\ \quad b &=& \frac{ \displaystyle \sum_{i=1}^{N} x_i \sum_{i=1}^{N} y_i - N \sum_{i=1}^{N} x_i y_i }{ \left( \displaystyle \sum_{i=1}^{N} x_i \right)^2 - N \displaystyle \sum_{i=1}^{N} x_i^2 } \end{eqnarray}\)

誤差最小値における\(a\)の導出(原点を通る制約のもと)

最小化する対象である誤差の2乗和は以下である。
\(\begin{eqnarray} \quad \displaystyle E &=& \sum_{i=1}^{N} (y_i - a x_i)^2 \\ \quad &=& \sum_{i=1}^{N} y_i^2-2 a x_i y_i+a^2 x_i^2 \\ \quad &=& \sum_{i=1}^{N} y_i^2-2 a \sum_{i=1}^{N} x_i y_i+a^2 \sum_{i=1}^{N} x_i^2 \end{eqnarray} \)
\(E\)が最小化であるとき、\(a\)について偏微分は0になる。
\( \quad \displaystyle \frac{\partial E}{\partial a} = -2 \sum_{i=1}^{N} x_i y_i + 2 a \sum_{i=1}^{N} x_i^2 = 0 \quad \rightarrow a = \displaystyle \sum_{i=1}^{N} x_i y_i / \displaystyle \sum_{i=1}^{N} x_i^2 \)

ソースコード


namespace DataFitting {

    /// <summary>線形フィッティング</summary>
    public class LinearFittingMethod : FittingMethod {

        /// <summary>コンストラクタ</summary>
        public LinearFittingMethod(FittingData[] data_list, bool is_enable_section) : base(data_list, is_enable_section ? 2 : 1) {
            IsEnableSection = is_enable_section;
        }

        /// <summary>y切片を有効にするか</summary>
        public bool IsEnableSection { get; private set; }

        /// <summary>フィッティング値</summary>
        public override double FittingValue(double x, Vector parameters) {
            if(parameters == null) {
                throw new ArgumentNullException(nameof(parameters));
            }
            if(parameters.Dim != ParametersCount) {
                throw new ArgumentException(nameof(parameters));
            }

            if(IsEnableSection) {
                return parameters[0] + parameters[1] * x;
            }
            else {
                return parameters[0] * x;
            }
        }

        /// <summary>フィッティング</summary>
        public Vector ExecureFitting() {
            if(IsEnableSection) {
                FittingData data;
                double sum_x = 0, sum_y = 0, sum_sq_x = 0, sum_xy = 0, n = data_list.Length;

                for(int i = 0; i < data_list.Length; i++) {
                    data = data_list[i];
                    sum_x += data.X;
                    sum_y += data.Y;
                    sum_sq_x += data.X * data.X;
                    sum_xy += data.X * data.Y;
                }

                double r = 1 / (sum_x * sum_x - n * sum_sq_x);
                double a = (sum_x * sum_xy - sum_sq_x * sum_y) * r;
                double b = (sum_x * sum_y - n * sum_xy) * r;

                return new Vector(a, b);
            }
            else {
                FittingData data;
                double sum_sq_x = 0, sum_xy = 0;

                for(int i = 0; i < data_list.Length; i++) {
                    data = data_list[i];
                    sum_sq_x += data.X * data.X;
                    sum_xy += data.X * data.Y;
                }

                return new Vector(sum_xy / sum_sq_x);
            }
        }

    }
}

単体テスト


namespace DataFitting.Tests {
    [TestClass()]
    public class LinearFittingMethodTests {
        [TestMethod()]
        public void ExecureFittingTest() {
            FittingData[] data_list = { new FittingData(2, 1), new FittingData(3, 8) };

            LinearFittingMethod fitting1 = new LinearFittingMethod(data_list, true);
            LinearFittingMethod fitting2 = new LinearFittingMethod(data_list, false);

            Assert.AreEqual(fitting1.ExecureFitting(), new Vector(-13, 7));
            Assert.AreEqual(fitting2.ExecureFitting(), new Vector(2));
        }
    }
}