関数fitting_line マニュアル

(The documentation of function fitting_line)
Last Update: 2021/12/7


◆機能・用途(Purpose)

データを直線でフィットする。
Fit a data by a straight line.

◆形式(Format)

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


◆引数(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.


◆戻り値(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 ∗line;
line=fitting_line(5,x,y);


◆計算式とアルゴリズム (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}\] と書ける。このとき \[\begin{equation} \myvector{G^T}\myvector{G}= \begin{pmatrix} \sum_{i=1}^N x_i^2 & \sum_{i=1}^N x_i \\ \sum_{i=1}^N x_i & N \end{pmatrix} \label{eq.GTG} \end{equation}\] \[\begin{equation} \myvector{G^T}\myvector{d}= \begin{pmatrix} \sum_{i=1}^N x_iy_i \\ \sum_{i=1}^N y_i \end{pmatrix} \label{eq.GTd} \end{equation}\] であり、最小二乗解は \[\begin{equation} \myvector{m}=(\myvector{G^T}\myvector{G})^{-1}\myvector{G^T}\myvector{d} \label{eq.solution} \end{equation}\] と書けるので、クラメルの公式を用いれば \[\begin{equation} a=\frac{ \begin{vmatrix} \sum_{i=1}^N x_iy_i & \sum_{i=1}^N x_i \\ \sum_{i=1}^N y_i & N \end{vmatrix} }{ \begin{vmatrix} \sum_{i=1}^N x_i^2 & \sum_{i=1}^N x_i \\ \sum_{i=1}^N x_i & N \end{vmatrix} } =\frac{1}{\Delta}\left[ N\sum_{i=1}^N x_iy_i -\left(\sum_{i=1}^N x_i\right)\left(\sum_{i=1}^N y_i\right) \right] \label{eq.a} \end{equation}\] \[\begin{equation} b=\frac{ \begin{vmatrix} \sum_{i=1}^N x_i^2 & \sum_{i=1}^N x_iy_i \\ \sum_{i=1}^N x_i & \sum_{i=1}^N y_i \end{vmatrix} }{ \begin{vmatrix} \sum_{i=1}^N x_i^2 & \sum_{i=1}^N x_i \\ \sum_{i=1}^N x_i & N \end{vmatrix} } =\frac{1}{\Delta}\left[ \left(\sum_{i=1}^N x_i^2\right)\left(\sum_{i=1}^N y_i\right) -\left(\sum_{i=1}^N x_i\right)\left(\sum_{i=1}^N x_iy_i\right) \right] \label{eq.b} \end{equation}\] を得る。ここで \[\begin{equation} \Delta = N\sum_{i=1}^N x_i^2-\left(\sum_{i=1}^N x_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}). Then Eqs. (\ref{eq.GTG}) and (\ref{eq.GTd}) hold, and the least squares solution 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 (\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[ N\sum_{i=1}^N x_i\delta_{ij} -\left(\sum_{i=1}^N x_i\right)\left(\sum_{i=1}^N \delta_{ij}\right) \right] \nonumber \\ &=& \frac{1}{\Delta}\left(Nx_j-\sum_{i=1}^N x_i\right) \label{eq.da.dyj} \end{eqnarray}\] \[\begin{eqnarray} \PartialDiff{b}{y_j} &=& \frac{1}{\Delta}\left[ \left(\sum_{i=1}^N x_i^2\right)\left(\sum_{i=1}^N \delta_{ij}\right) -\left(\sum_{i=1}^N x_i\right)\left(\sum_{i=1}^N x_i\delta_{ij}\right) \right] \nonumber \\ &=& \frac{1}{\Delta}\left(\sum_{i=1}^N x_i^2-x_j\sum_{i=1}^N x_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(Nx_j-\sum_{i=1}^N x_i\right)^2e_j^2} \nonumber \\ &=& \frac{1}{\Delta}\sqrt{\sum_{j=1}^N \left(Nx_j-\sum_{i=1}^N x_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(\sum_{i=1}^N x_i^2-x_j\sum_{i=1}^N x_i\right)^2e_j^2} \nonumber \\ &=& \frac{1}{\Delta}\sqrt{\sum_{j=1}^N \left(\sum_{i=1}^N x_i^2-x_j\sum_{i=1}^N x_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 given by 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.