重み付き多項式フィッティング

最適化
重み付き多項式フィッティング; Weighted Polynomial Fitting; Robust Estimation

概要

多項式フィッティングとは独立変数\(x_i\)と従属変数\(y_i\)とフィッティングに対する重み\(w_i\)の組で表される\(N\)個のデータ\(i = 1 \cdots N \)があるとき、
\(n\)次多項式\(\tilde{y_i} = \displaystyle \sum_{k=0}^{n} r_k x_i^k\)と\(y_i\)との誤差と重み\(w_i\)の積和を最小化する、\(\boldsymbol{r} = (r_0 \cdots r_n)^T \)を求めることである。
誤差には2乗和が用いられる。
\(\quad \displaystyle E = \sum_{i=1}^{N} w_i (y_i - \tilde{y_i})^2 \)

なお、データの誤差の分散\(\sigma^2 \)が既知の場合、重み\(w\)は以下で与えられる。
\(\quad \displaystyle w = \frac{1}{{\sigma}^2} \)

重み付き多項式フィッティング

誤差最小値における\(\boldsymbol{r} \)の導出

最小化する対象である誤差の2乗和は以下である。
\(\begin{eqnarray} \quad \displaystyle E &=& \sum_{i=1}^{N} w_i (y_i - \sum_{k=0}^{n} r_k x_i^k)^2 \\ \quad &=& | W (\boldsymbol{y} - M \boldsymbol{r})^2 | \quad \left( \boldsymbol{y} = (y_1 \cdots y_N)^T, \ \ W_{ij}= \begin{cases} w_i &\quad i=j \\ 0 &\quad i \neq j \end{cases}, \ \ M_{ki} = x_i^k \right) \\ \quad &=& W ( \boldsymbol{y} - M \boldsymbol{r})^T (\boldsymbol{y} - M \boldsymbol{r}) & \cdots E \geq 0 \\ \quad &=& ( \boldsymbol{y}^T - \boldsymbol{r}^T M^T) (W \boldsymbol{y} - W M \boldsymbol{r}) \\ \quad &=& \boldsymbol{y}^T W \boldsymbol{y} -\boldsymbol{y}^T W M \boldsymbol{r} - \boldsymbol{r}^T M^T W \boldsymbol{y} + \boldsymbol{r}^T M^T W M \boldsymbol{r} \\ \quad &=& \boldsymbol{y}^T W \boldsymbol{y} - (\boldsymbol{y}^T W M \boldsymbol{r})^T - \boldsymbol{r}^T M^T W \boldsymbol{y} + \boldsymbol{r}^T M^T W M \boldsymbol{r} & \cdots \boldsymbol{y}^T W M \boldsymbol{r} \in \mathbb{R}\\ \quad &=& \boldsymbol{y}^T W \boldsymbol{y} - 2 \boldsymbol{r}^T M^T W \boldsymbol{y} + \boldsymbol{r}^T M^T W M \boldsymbol{r} & \cdots W^T=W \\ \end{eqnarray} \)
\(E\)が最小化であるとき、\(\boldsymbol{r}\)についての偏微分は\(\boldsymbol{0}\)になる。
\( \quad \displaystyle \frac{\partial E}{\partial \boldsymbol{r}} = - 2 M^T W \boldsymbol{y} + 2 M^T W M \boldsymbol{r} = \boldsymbol{0} \\ \quad \rightarrow M^T W M \boldsymbol{r} = M^T W \boldsymbol{y} \\ \quad \rightarrow \boldsymbol{r} = (M^T W M)^{-1} M^T W \boldsymbol{y} \\ \)
\((M^T W M)^{-1}\)ではパラメータの数(\(n+1\))の正方行列の逆行列を計算することになる。
原点を通る多項式を得たい場合は\(\displaystyle \sum_{k=1}^{n} r_k x_i^k\)としてパラメータの数を1減らし、その後の導出は同様になる。

ソースコード


namespace DataFitting {

    /// <summary>重み付き多項式フィッティング</summary>
    public class WeightedPolynomialFittingMethod : FittingMethod {

        readonly Matrix weight_matrix;

        /// <summary>コンストラクタ</summary>
        public WeightedPolynomialFittingMethod(FittingData[] data_list, double[] weight_list, int degree, bool is_enable_section)
            : base(data_list, degree + (is_enable_section ? 1 : 0)) {

            this.Degree = degree;
            this.IsEnableSection = is_enable_section;

            if(weight_list == null) {
                throw new ArgumentNullException(nameof(weight_list));
            }

            if(data_list.Length != weight_list.Length) {
                throw new ArgumentException($"{nameof(data_list)},{nameof(weight_list)}");
            }

            foreach(var weight in weight_list) {
                if(!(weight >= 0)) {
                    throw new ArgumentException(nameof(weight_list));
                }
            }

            this.weight_matrix = new Matrix(weight_list.Length, weight_list.Length);

            for(int i = 0; i < weight_list.Length; i++) {
                weight_matrix[i, i] = weight_list[i];
            }
        }

        /// <summary>次数</summary>
        public int Degree {
            get; private set;
        }

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

        /// <summary>重み付き誤差二乗和</summary>
        public double WeightedCost(Vector coefficients) {
            if(coefficients == null) {
                throw new ArgumentNullException(nameof(coefficients));
            }
            if(coefficients.Dim != ParametersCount) {
                throw new ArgumentException(nameof(coefficients));
            }

            Vector errors = Error(coefficients);
            double cost = 0;
            for(int i = 0; i < errors.Dim; i++) {
                cost += weight_matrix[i, i] * errors[i] * errors[i];
            }

            return cost;
        }

        /// <summary>フィッティング値</summary>
        public override double FittingValue(double x, Vector coefficients) {
            if(IsEnableSection) {
                double y = coefficients[0], ploy_x = 1;

                for(int i = 1; i < coefficients.Dim; i++) {
                    ploy_x *= x;
                    y += ploy_x * coefficients[i];
                }

                return y;
            }
            else {
                double y = 0, ploy_x = 1;

                for(int i = 0; i < coefficients.Dim; i++) {
                    ploy_x *= x;
                    y += ploy_x * coefficients[i];
                }

                return y;
            }
        }

        /// <summary>フィッティング</summary>
        public Vector ExecureFitting() {
            Matrix m = new Matrix(data_list.Length, ParametersCount);
            Vector b = Vector.Zero(data_list.Length);

            if(IsEnableSection) {
                for(int i = 0; i < data_list.Length; i++) {
                    double x = data_list[i].X;
                    b[i] = data_list[i].Y;

                    m[i, 0] = 1;

                    for(int j = 1; j <= Degree; j++) {
                        m[i, j] = m[i, j - 1] * x;
                    }
                }
            }
            else {
                for(int i = 0; i < data_list.Length; i++) {
                    double x = data_list[i].X;
                    b[i] = data_list[i].Y;

                    m[i, 0] = x;

                    for(int j = 1; j < Degree; j++) {
                        m[i, j] = m[i, j - 1] * x;
                    }
                }
            }

            return (m.Transpose * weight_matrix * m).Inverse * m.Transpose * weight_matrix * b;
        }
    }
}