最適化重み付き直線フィッティング; 重み付き線形回帰; Weighted Linear Regression
概要
直線フィッティングとは独立変数\(x_i\)と従属変数\(y_i\)とフィッティングに対する重み\(w_i\)の組で表される\(N\)個のデータ\(i = 1 \cdots N \)があるとき、
直線\(\tilde{y_i} = a + b x_i\)または原点を通る直線\(\tilde{y_i} = a x_i\)と、\(y_i\)との誤差と重み\(w_i\)の積和を最小化する、\(a, b\)の組または\(a\)を求めることである。
誤差には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} \)
誤差最小値における\(a, b\)の導出
最小化する対象である誤差の2乗和は以下である。
\(\begin{eqnarray} \quad \displaystyle E &=& \sum_{i=1}^{N} w_i (y_i - a - b x_i)^2 \\ \quad &=& \sum_{i=1}^{N} w_i y_i^2-2 b w_i x_i y_i-2 a w_i y_i+b^2 w_i x_i^2+2 a b w_i x_i+a^2 w_i \\ \quad &=& \sum_{i=1}^{N} w_i y_i^2-2 b \sum_{i=1}^{N} w_i x_i y_i-2 a \sum_{i=1}^{N} w_i y_i+b^2 \sum_{i=1}^{N} w_i x_i^2+2 a b \sum_{i=1}^{N} w_i x_i+ a^2 \sum_{i=1}^{N} w_i \\ \end{eqnarray} \)
\(E\)が最小化であるとき、\(a, b\)について偏微分はともに0になる。
\(\begin{eqnarray} \quad \displaystyle \frac{\partial E}{\partial a} &=& -2 \sum_{i=1}^{N} w_i y_i+2 b \sum_{i=1}^{N} w_i x_i+2 a \sum_{i=1}^{N} w_i = 0 \quad &\rightarrow& a \sum_{i=1}^{N} w_i + b \sum_{i=1}^{N} w_i x_i = \sum_{i=1}^{N} w_i y_i \\ \quad \displaystyle \frac{\partial E}{\partial b} &=& -2 \sum_{i=1}^{N} w_i x_i y_i + 2 b \sum_{i=1}^{N} w_i x_i^2 + 2 a \sum_{i=1}^{N} w_i x_i = 0 \quad &\rightarrow& a \sum_{i=1}^{N} w_i x_i + b \sum_{i=1}^{N} w_i x_i^2 = \sum_{i=1}^{N} w_i x_i y_i \\ \end{eqnarray}\)
2式を連立させて\(a, b\)について解く。
\(\quad \begin{pmatrix} \sum_{i=1}^{N} w_i & \sum_{i=1}^{N} w_i x_i \\ \sum_{i=1}^{N} w_i x_i & \sum_{i=1}^{N} w_i x_i^2 \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} \sum_{i=1}^{N} w_i y_i \\ \sum_{i=1}^{N} w_i x_i y_i \end{pmatrix} \)
であるから
\( \begin{eqnarray} \quad a &=& \frac{ \displaystyle \sum_{i=1}^{N} w_i x_i \sum_{i=1}^{N} w_i x_i y_i - \sum_{i=1}^{N} w_i x_i^2 \sum_{i=1}^{N} w_i y_i }{ \left( \displaystyle \sum_{i=1}^{N} w_i x_i \right)^2 - \displaystyle \sum_{i=1}^{N} w_i \displaystyle \sum_{i=1}^{N} w_i x_i^2 } \\ \quad b &=& \frac{ \displaystyle \sum_{i=1}^{N} w_i x_i \sum_{i=1}^{N} w_i y_i - \sum_{i=1}^{N} w_i \sum_{i=1}^{N} w_i x_i y_i }{ \left( \displaystyle \sum_{i=1}^{N} w_i x_i \right)^2 - \displaystyle \sum_{i=1}^{N} w_i \displaystyle \sum_{i=1}^{N} w_i x_i^2 } \end{eqnarray}\)
誤差最小値における\(a\)の導出(原点を通る制約のもと)
最小化する対象である誤差の2乗和は以下である。
\(\begin{eqnarray} \quad \displaystyle E &=& \sum_{i=1}^{N} w_i (y_i - a x_i)^2 \\ \quad &=& \sum_{i=1}^{N} w_i y_i^2-2 a w_i x_i y_i+a^2 w_i x_i^2 \\ \quad &=& \sum_{i=1}^{N} w_i y_i^2-2 a \sum_{i=1}^{N} w_i x_i y_i+a^2 \sum_{i=1}^{N} w_i x_i^2 \end{eqnarray} \)
\(E\)が最小化であるとき、\(a\)について偏微分は0になる。
\( \quad \displaystyle \frac{\partial E}{\partial a} = -2 \sum_{i=1}^{N} w_i x_i y_i + 2 a \sum_{i=1}^{N} w_i x_i^2 = 0 \quad \rightarrow a = \displaystyle \sum_{i=1}^{N} w_i x_i y_i / \displaystyle \sum_{i=1}^{N} w_i x_i^2 \)
ソースコード
namespace DataFitting {
/// <summary>重み付き線形フィッティング</summary>
public class WeightedLinearFittingMethod : FittingMethod {
readonly double[] weight_list;
/// <summary>コンストラクタ</summary>
public WeightedLinearFittingMethod(FittingData[] data_list, double[] weight_list, bool is_enable_section) : base(data_list, is_enable_section ? 2 : 1) {
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_list = weight_list;
}
/// <summary>y切片を有効にするか</summary>
public bool IsEnableSection { get; private set; }
/// <summary>重み付き誤差二乗和</summary>
public double WeightedCost(Vector parameters) {
if(parameters == null) {
throw new ArgumentNullException(nameof(parameters));
}
if(parameters.Dim != ParametersCount) {
throw new ArgumentException(nameof(parameters));
}
Vector errors = Error(parameters);
double cost = 0;
for(int i = 0; i < errors.Dim; i++) {
cost += weight_list[i] * errors[i] * errors[i];
}
return cost;
}
/// <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(weight_list == null) {
throw new InvalidOperationException();
}
if(IsEnableSection) {
FittingData data;
double w, sum_w = 0, sum_wx = 0, sum_wy = 0, sum_wxx = 0, sum_wxy = 0;
for(int i = 0; i < data_list.Length; i++) {
data = data_list[i];
w = weight_list[i];
sum_w += w;
sum_wx += w * data.X;
sum_wy += w * data.Y;
sum_wxx += w * data.X * data.X;
sum_wxy += w * data.X * data.Y;
}
double r = 1 / (sum_wx * sum_wx - sum_w * sum_wxx);
double a = (sum_wx * sum_wxy - sum_wxx * sum_wy) * r;
double b = (sum_wx * sum_wy - sum_w * sum_wxy) * r;
return new Vector(a, b);
}
else {
FittingData data;
double w, sum_wxx = 0, sum_wxy = 0;
for(int i = 0; i < data_list.Length; i++) {
data = data_list[i];
w = weight_list[i];
sum_wxx += w * data.X * data.X;
sum_wxy += w * data.X * data.Y;
}
return new Vector(sum_wxy / sum_wxx);
}
}
}
}
関連項目
関数フィッティング手法基本クラス
ベクトルクラス
行列クラス
直線フィッティング
多項式フィッティング
非線形フィッティング Gauss-Newton法
非線形フィッティング Levenberg-Marquardt法
重み付き多項式フィッティング
最小二乗法におけるロバスト推定
