関数fitting_line_weight マニュアル

(The documentation of function fitting_line_weight)
Last Update: 2021/12/1


◆機能・用途(Purpose)

データを直線で重み付きフィットする。
Fit a data by a straight line with weights.

◆形式(Format)

#include <fitting.h>
inline double ∗fitting_line_weight
(const int N,const double ∗x,const double ∗y, const double ∗weight)


◆引数(Arguments)

N フィットするデータ点数。
The number of data points to fit.
x フィットするデータの\(x\)の値のリスト。 \(N\)個の要素から成る配列で与える。
The list of data values of \(x\) to fit, given as an array composed of \(N\) components.
y フィットするデータの\(y\)の値のリスト。 \(N\)個の要素から成る配列で与える。
The list of data values of \(y\) to fit, given as an array composed of \(N\) components.
weight 各データに対する重み係数のリスト。 \(N\)個の要素から成る配列で与える。
The list of weights for individual data, given as an array of \(N\) components.


◆戻り値(Return value)

直線の傾きを第1要素、切片を第2要素、 傾きの推定誤差を第3要素、切片の推定誤差を第4要素とする配列。
An array composed of the slope, intercept, and estimation errors of them as the 1st to 4th components, respectively.

◆使用例(Example)

double x[]={1.0,2.0,3.0,4.0,5.0};
double y[]={1.1,2.1,3.4,4.5,5.5};
double weight[]={1.0,0.5,1.0,2.0,1.0,1.5};
double ∗line;
line=fitting_line_weight(5,x,y,weight);


◆計算式とアルゴリズム (Formula and algorithm)

1. 傾きと切片の式 (Equations for the slope and intercept)

\((x,y)\)平面上の\(N\)個のデータ点\((x_1,y_1),\cdots,(x_N,y_N)\)を 直線\(y=ax+b\)でフィットすることを考える。 観測方程式は \[\begin{equation} y_i=ax_i+b \hspace{1em} (i=1,\cdots,N) \label{eq.obs} \end{equation}\] であるので行列形式で \[\begin{equation} \myvector{d}=\myvector{G}\myvector{m} \label{eq.obs.matrix} \end{equation}\] \[\begin{equation} \myvector{d}= \begin{pmatrix} y_1 \\ \vdots \\ y_N \end{pmatrix} \label{eq.d} \end{equation}\] \[\begin{equation} \myvector{G}= \begin{pmatrix} x_1 & 1 \\ \vdots & \vdots \\ x_N & 1 \end{pmatrix} \label{eq.G} \end{equation}\] \[\begin{equation} \myvector{m}= \begin{pmatrix} a \\ b \end{pmatrix} \label{eq.m} \end{equation}\] と書ける。 各データに対する重み係数を\(w_1,\cdots,w_N\)として \[\begin{equation} \myvector{W}= \begin{pmatrix} w_1 & 0 & \cdots & 0 \\ 0 & w_2 & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & w_N \end{pmatrix} \label{eq.W} \end{equation}\] とおく。このとき \[\begin{equation} \myvector{G^T}\myvector{W}\myvector{G}= \begin{pmatrix} \sum_{i=1}^N w_ix_i^2 & \sum_{i=1}^N w_ix_i \\ \sum_{i=1}^N w_ix_i & \sum_{i=1}^N w_i \end{pmatrix} \label{eq.GTG} \end{equation}\] \[\begin{equation} \myvector{G^T}\myvector{W}\myvector{d}= \begin{pmatrix} \sum_{i=1}^N w_ix_iy_i \\ \sum_{i=1}^N w_iy_i \end{pmatrix} \label{eq.GTd} \end{equation}\] であり、 \[\begin{equation} E= \left(\myvector{d}-\myvector{G}\myvector{m}\right)^T \myvector{W} \left(\myvector{d}-\myvector{G}\myvector{m}\right) \label{eq.E} \end{equation}\] を最小にする重み付き最小二乗解は \[\begin{equation} \myvector{m} =(\myvector{G^T}\myvector{W}\myvector{G})^{-1} \myvector{G^T}\myvector{W}\myvector{d} \label{eq.solution} \end{equation}\] と書けるので、クラメルの公式を用いれば \[\begin{eqnarray} a &=& \frac{ \begin{vmatrix} \sum_{i=1}^N w_ix_iy_i & \sum_{i=1}^N w_ix_i \\ \sum_{i=1}^N w_iy_i & \sum_{i=1}^N w_i \end{vmatrix} }{ \begin{vmatrix} \sum_{i=1}^N w_ix_i^2 & \sum_{i=1}^N w_ix_i \\ \sum_{i=1}^N w_ix_i & \sum_{i=1}^N w_i \end{vmatrix} } \nonumber \\ &=& \frac{1}{\Delta}\left[ \left(\sum_{i=1}^N w_ix_iy_i\right)\left(\sum_{i=1}^N w_i\right) -\left(\sum_{i=1}^N w_ix_i\right)\left(\sum_{i=1}^N w_iy_i\right) \right] \label{eq.a} \end{eqnarray}\] \[\begin{eqnarray} b &=& \frac{ \begin{vmatrix} \sum_{i=1}^N w_ix_i^2 & \sum_{i=1}^N w_ix_iy_i \\ \sum_{i=1}^N w_ix_i & \sum_{i=1}^N w_iy_i \end{vmatrix} }{ \begin{vmatrix} \sum_{i=1}^N w_ix_i^2 & \sum_{i=1}^N w_ix_i \\ \sum_{i=1}^N w_ix_i & \sum_{i=1}^N w_i \end{vmatrix} } \nonumber \\ &=& \frac{1}{\Delta}\left[ \left(\sum_{i=1}^N w_ix_i^2\right)\left(\sum_{i=1}^N w_iy_i\right) -\left(\sum_{i=1}^N w_ix_i\right)\left(\sum_{i=1}^N w_ix_iy_i\right) \right] \label{eq.b} \end{eqnarray}\] を得る。ここで \[\begin{equation} \Delta= \left(\sum_{i=1}^N w_ix_i^2\right)\left(\sum_{i=1}^N w_i\right) -\left(\sum_{i=1}^N w_ix_i\right)^2 \label{eq.Delta} \end{equation}\] とおいた。
Consider a problem to fit \(N\) data points on \((x,y)\) plane, \((x_1,y_1),\cdots,(x_N,y_N)\), by a straight line \(y=ax+b\). The observation equation is given by Eq. (\ref{eq.obs}), which can be arranged in matrix form as Eqs. (\ref{eq.obs.matrix})-(\ref{eq.m}). Let \(w_1,\cdots,w_N\) be the weight for each data and a matrix \(W\) be defined by Eq. (\ref{eq.W}). Then Eqs. (\ref{eq.GTG}) and (\ref{eq.GTd}) hold, and the weighted least squares solution, which minimizes \(E\) defined by Eq. (\ref{eq.E}), is given by Eq. (\ref{eq.solution}). Using Cramer's formula, one obtains Eqs. (\ref{eq.a}) and (\ref{eq.b}), where \(\Delta\) is defined by eq. (\ref{eq.Delta}).

2. 推定誤差の式 (Equations for the estimation errors)

\(x_i\)は誤差を持たないものとし、\(y_i\)の誤差はデータと直線との残差 \[\begin{equation} e_i = y_i-(ax_i+b) = y_i-ax_i-b \label{eq.ei} \end{equation}\] で近似する。このとき誤差伝搬則より\(a\), \(b\)の推定誤差はそれぞれ \[\begin{equation} \sigma_a=\sqrt{\sum_{j=1}^N \left(\PartialDiff{a}{y_j}\right)^2e_j^2} \label{eq.sigma_a.original} \end{equation}\] \[\begin{equation} \sigma_b=\sqrt{\sum_{j=1}^N \left(\PartialDiff{b}{y_j}\right)^2e_j^2} \label{eq.sigma_b.original} \end{equation}\] と書ける。(\ref{eq.a})(\ref{eq.b})式および \(\PartialDiff{y_i}{y_j}=\delta_{ij}\)より \[\begin{eqnarray} \PartialDiff{a}{y_j} &=& \frac{1}{\Delta}\left[ \left(\sum_{i=1}^N w_ix_i\delta_{ij}\right)\left(\sum_{i=1}^N w_i\right) -\left(\sum_{i=1}^N w_ix_i\right)\left(\sum_{i=1}^N w_i\delta_{ij}\right) \right] \nonumber \\ &=& \frac{1}{\Delta} \left(w_jx_j\sum_{i=1}^N w_i-w_j\sum_{i=1}^N w_ix_i\right) \label{eq.da.dyj} \end{eqnarray}\] \[\begin{eqnarray} \PartialDiff{b}{y_j} &=& \frac{1}{\Delta}\left[ \left(\sum_{i=1}^N w_ix_i^2\right) \left(\sum_{i=1}^N w_i\delta_{ij}\right) -\left(\sum_{i=1}^N w_ix_i\right) \left(\sum_{i=1}^N w_ix_i\delta_{ij}\right) \right] \nonumber \\ &=& \frac{1}{\Delta} \left(w_j\sum_{i=1}^N w_ix_i^2-w_jx_j\sum_{i=1}^N w_ix_i\right) \label{eq.db.dyj} \end{eqnarray}\] であるので、これらを (\ref{eq.sigma_a.original})(\ref{eq.sigma_b.original})に代入して \[\begin{eqnarray} \sigma_a &=& \sqrt{\sum_{j=1}^N \frac{1}{\Delta^2} \left(w_jx_j\sum_{i=1}^N w_i-w_j\sum_{i=1}^N w_ix_i\right)^2 e_j^2} \nonumber \\ &=& \frac{1}{\Delta}\sqrt{\sum_{j=1}^N \left(w_jx_j\sum_{i=1}^N w_i-w_j\sum_{i=1}^N w_ix_i\right)^2 \left(y_j-ax_j-b\right)^2} \label{eq.sigma_a} \end{eqnarray}\] \[\begin{eqnarray} \sigma_b &=& \sqrt{\sum_{j=1}^N \frac{1}{\Delta^2} \left(w_j\sum_{i=1}^N w_ix_i^2-w_jx_j\sum_{i=1}^N w_ix_i\right)^2 e_j^2} \nonumber \\ &=& \frac{1}{\Delta}\sqrt{\sum_{j=1}^N \left(w_j\sum_{i=1}^N w_ix_i^2-w_jx_j\sum_{i=1}^N w_ix_i\right)^2 \left(y_j-ax_j-b\right)^2} \label{eq.sigma_b} \end{eqnarray}\] を得る。
No error is assumed for \(x_i\), and the error of \(y_i\) is approximated by the residual between the data and the straight line (eq. \ref{eq.ei}). Then the error propagation law tells that the estimation errors for \(a\) and \(b\) are eqs (\ref{eq.sigma_a.original}) and \(\ref{eq.sigma_b.original}\), respectively. From eqs (\ref{eq.a}), (\ref{eq.b}), and a relation \(\PartialDiff{y_i}{y_j}=\delta_{ij}\), one obtains (\ref{eq.da.dyj}) and (\ref{eq.db.dyj}). Inserting these relations into (\ref{eq.sigma_a.original}) and \(\ref{eq.sigma_b.original}\) results in (\ref{eq.sigma_a}) and \(\ref{eq.sigma_b}\), respectively.