関数fitting_quadratic マニュアル

(The documentation of function fitting_quadratic)
Last Update: 2022/12/19


◆機能・用途(Purpose)

データを2次関数でフィットする。
Fit a data by a quadractic function.

◆形式(Format)

#include <fitting.h>
inline double ∗fitting_quadratic
(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)

データを2次関数\(y=ax^2+bx+c\)でフィットしたときの \(a\)を第1要素、\(b\)を第2要素、\(c\)を第3要素、 それらの推定誤差を第4-第6要素とする配列。
An array composed of the estimated values of \(a\), \(b\), \(c\) in the 1st to 3rd components and their estimation errors in the 4th to 6th components.

◆使用例(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 ∗quadratic;
quadratic=fitting_quadratic(5,x,y);


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

1. \(a\), \(b\), \(c\)の最適値 (Optimal values of \(a\), \(b\), and \(c\))

\((x,y)\)平面上の\(N\)個のデータ点\((x_1,y_1),\cdots,(x_N,y_N)\)を 2次関数\(y=ax^2+bx+c\)でフィットすることを考える。 観測方程式は \[\begin{equation} y_i=ax_i^2+bx_i+c \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^2 & x_1 & 1 \\ \vdots & \vdots \\ x_N^2 & x_N & 1 \end{pmatrix} \label{eq.G} \end{equation}\] \[\begin{equation} \myvector{m}= \begin{pmatrix} a \\ b \\ c \end{pmatrix} \label{eq.m} \end{equation}\] と書ける。このとき \[\begin{equation} \myvector{G^T}\myvector{G}= \begin{pmatrix} \sum_{i=1}^N x_i^4 & \sum_{i=1}^N x_i^3 & \sum_{i=1}^N x_i^2 \\ \sum_{i=1}^N x_i^3 & \sum_{i=1}^N x_i^2 & \sum_{i=1}^N x_i \\ \sum_{i=1}^N x_i^2 & \sum_{i=1}^N x_i & N \end{pmatrix} =N \begin{pmatrix} <x^4> & <x^3> & <x^2> \\ <x^3> & <x^2> & <x^1> \\ <x^2> & <x^1> & <x^0> \end{pmatrix} \label{eq.GTG} \end{equation}\] \[\begin{equation} \myvector{G^T}\myvector{d}= \begin{pmatrix} \sum_{i=1}^N x_i^2y_i \\ \sum_{i=1}^N x_iy_i \\ \sum_{i=1}^N y_i \end{pmatrix} =N \begin{pmatrix} <x^2y> \\ <xy> \\ <y> \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{eqnarray} a &=& \frac{ \begin{vmatrix} <x^2y> & <x^3> & <x^2> \\ <xy> & <x^2> & <x^1> \\ <y> & <x^1> & <x^0> \end{vmatrix} }{ \begin{vmatrix} <x^4> & <x^3> & <x^2> \\ <x^3> & <x^2> & <x^1> \\ <x^2> & <x^1> & <x^0> \end{vmatrix} } \nonumber \\ &=& \frac{1}{\Delta} \left[ \left( <x^2> <x^0> - <x^1>^2 \right) <x^2y> \right. \nonumber \\ & & -\left( <x^3> <x^0> - <x^2> <x^1> \right) <xy> \nonumber \\ & & \left. +\left( <x^3> <x^1> - <x^2>^2 \right) <y> \right] \label{eq.a} \end{eqnarray}\] \[\begin{eqnarray} b &=& \frac{ \begin{vmatrix} <x^4> & <x^2y> & <x^2> \\ <x^3> & <xy> & <x^1> \\ <x^2> & <y> & <x^0> \end{vmatrix} }{ \begin{vmatrix} <x^4> & <x^3> & <x^2> \\ <x^3> & <x^2> & <x^1> \\ <x^2> & <x^1> & <x^0> \end{vmatrix} } \nonumber \\ &=& \frac{1}{\Delta} \left[ -\left( <x^3> <x^0> - <x^2> <x^1> \right) <x^2y> \right. \nonumber \\ & & +\left( <x^4> <x^0> - <x^2>^2 \right) <xy> \nonumber \\ & & \left. -\left( <x^4> <x^1> - <x^3> <x^2> \right) <y> \right] \label{eq.b} \end{eqnarray}\] \[\begin{eqnarray} c &=& \frac{ \begin{vmatrix} <x^4> & <x^3> & <x^2y> \\ <x^3> & <x^2> & <xy> \\ <x^2> & <x^1> & <y> \end{vmatrix} }{ \begin{vmatrix} <x^4> & <x^3> & <x^2> \\ <x^3> & <x^2> & <x^1> \\ <x^2> & <x^1> & <x^0> \end{vmatrix} } \nonumber \\ &=& \frac{1}{\Delta} \left[ \left( <x^3> <x^1> - <x^2>^2 \right) <x^2y> \right. \nonumber \\ & & -\left( <x^4> <x^1> - <x^3> <x^2> \right) <xy> \nonumber \\ & & \left. +\left( <x^4> <x^2> - <x^3>^2 \right) <y> \right] \label{eq.c} \end{eqnarray}\] を得る。ここで \[\begin{eqnarray} \Delta &=& <x^4> <x^2> <x^0> +2 <x^3> <x^2> <x^1> \nonumber \\ & & - <x^4> <x^1>^2 - <x^3>^2 <x^0> - <x^2>^3 \label{eq.Delta} \end{eqnarray}\] とおいた。
Consider a problem to fit \(N\) data points on \((x,y)\) plane, \((x_1,y_1),\cdots,(x_N,y_N)\), by a quadratic function \(y=ax^2+bx+c\). 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, where < > represents an average, and the least squares solution is given by Eq. (\ref{eq.solution}). Using Cramer's formula, one obtains Eqs. (\ref{eq.a})-(\ref{eq.c}), where \(\Delta\) is defined by (\ref{eq.Delta}).

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

\(x_i\)は誤差を持たないものとし、\(y_i\)の誤差を\(\sigma_{yi}\)とする。 このとき誤差伝搬則より\(a\), \(b\), \(c\)の推定誤差はそれぞれ \[\begin{equation} \sigma_a=\sqrt{\sum_{j=1}^N \left(\PartialDiff{a}{y_j}\right)^2\sigma_{yj}^2} \label{eq.sigma_a.original} \end{equation}\] \[\begin{equation} \sigma_b=\sqrt{\sum_{j=1}^N \left(\PartialDiff{b}{y_j}\right)^2\sigma_{yj}^2} \label{eq.sigma_b.original} \end{equation}\] \[\begin{equation} \sigma_c=\sqrt{\sum_{j=1}^N \left(\PartialDiff{c}{y_j}\right)^2\sigma_{yj}^2} \label{eq.sigma_c.original} \end{equation}\] と書ける。 (\ref{eq.a})-(\ref{eq.c})式より \[\begin{eqnarray} \PartialDiff{a}{y_j} &=& \frac{1}{\Delta} \left[ \left( <x^2> <x^0> - <x^1>^2 \right) \PartialDiff{<x^2y>}{y_j} \right. \nonumber \\ & & -\left( <x^3> <x^0> - <x^2> <x^1> \right) \PartialDiff{<xy>}{y_j} \nonumber \\ & & \left. +\left( <x^3> <x^1> - <x^2>^2 \right) \PartialDiff{<y>}{y_j} \right] \label{eq.da_dy.derive1} \end{eqnarray}\] \[\begin{eqnarray} \PartialDiff{b}{y_j} &=& \frac{1}{\Delta} \left[ -\left( <x^3> <x^0> - <x^2> <x^1> \right) \PartialDiff{<x^2y>}{y_j} \right. \nonumber \\ & & +\left( <x^4> <x^0> - <x^2>^2 \right) \PartialDiff{<xy>}{y_j} \nonumber \\ & & \left. -\left( <x^4> <x^1> - <x^3> <x^2> \right) \PartialDiff{<y>}{y_j} \right] \label{eq.db_dy.derive1} \end{eqnarray}\] \[\begin{eqnarray} \PartialDiff{c}{y_j} &=& \frac{1}{\Delta} \left[ \left( <x^3> <x^1> - <x^2>^2 \right) \PartialDiff{<x^2y>}{y_j} \right. \nonumber \\ & & -\left( <x^4> <x^1> - <x^3> <x^2> \right) \PartialDiff{<xy>}{y_j} \nonumber \\ & & \left. +\left( <x^4> <x^2> - <x^3>^2 \right) \PartialDiff{<y>}{y_j} \right] \label{eq.dc_dy.derive1} \end{eqnarray}\] である。 \[\begin{eqnarray} \PartialDiff{<x^2y>}{y_j} &=& \PartialDiff{}{y_j} \left(\frac{1}{N}\sum_{i=1}^N x_i^2y_i\right) \nonumber \\ &=& \frac{1}{N}\sum_{i=1}^N x_i^2\PartialDiff{y_i}{y_j} \nonumber \\ &=& \frac{1}{N}\sum_{i=1}^N x_i^2\delta_{ij} \nonumber \\ &=& \frac{1}{N}x_j^2 \label{eq.dx2y_dy} \end{eqnarray}\] \[\begin{eqnarray} \PartialDiff{<xy>}{y_j} &=& \PartialDiff{}{y_j} \left(\frac{1}{N}\sum_{i=1}^N x_iy_i\right) \nonumber \\ &=& \frac{1}{N}\sum_{i=1}^N x_i\PartialDiff{y_i}{y_j} \nonumber \\ &=& \frac{1}{N}\sum_{i=1}^N x_i\delta_{ij} \nonumber \\ &=& \frac{1}{N}x_j \label{eq.dxy_dy} \end{eqnarray}\] \[\begin{eqnarray} \PartialDiff{<y>}{y_j} &=& \PartialDiff{}{y_j} \left(\frac{1}{N}\sum_{i=1}^N y_i\right) \nonumber \\ &=& \frac{1}{N}\sum_{i=1}^N \PartialDiff{y_i}{y_j} \nonumber \\ &=& \frac{1}{N}\sum_{i=1}^N \delta_{ij} \nonumber \\ &=& \frac{1}{N} \label{eq.dy_dy} \end{eqnarray}\] であるので、 これらを(\ref{eq.da_dy.derive1})-(\ref{eq.dc_dy.derive1})に代入して \[\begin{eqnarray} \PartialDiff{a}{y_j} &=& \frac{1}{\Delta} \left[ \left( <x^2> <x^0> - <x^1>^2 \right) \frac{1}{N}x_j^2 \right. \nonumber \\ & & -\left( <x^3> <x^0> - <x^2> <x^1> \right) \frac{1}{N}x_j \nonumber \\ & & \left. +\left( <x^3> <x^1> - <x^2>^2 \right) \frac{1}{N} \right] \nonumber \\ &=& \frac{1}{N\Delta} \left[ \left( <x^2> <x^0> - <x^1>^2 \right) x_j^2 \right. \nonumber \\ & & -\left( <x^3> <x^0> - <x^2> <x^1> \right) x_j \nonumber \\ & & \left. +\left( <x^3> <x^1> - <x^2>^2 \right) \right] \label{eq.da_dy.derive2} \end{eqnarray}\] \[\begin{eqnarray} \PartialDiff{b}{y_j} &=& \frac{1}{\Delta} \left[ -\left( <x^3> <x^0> - <x^2> <x^1> \right) \frac{1}{N}x_j^2 \right. \nonumber \\ & & +\left( <x^4> <x^0> - <x^2>^2 \right) \frac{1}{N}x_j \nonumber \\ & & \left. -\left( <x^4> <x^1> - <x^3> <x^2> \right) \frac{1}{N} \right] \nonumber \\ &=& \frac{1}{N\Delta} \left[ -\left( <x^3> <x^0> - <x^2> <x^1> \right) x_j^2 \right. \nonumber \\ & & +\left( <x^4> <x^0> - <x^2>^2 \right) x_j \nonumber \\ & & \left. -\left( <x^4> <x^1> - <x^3> <x^2> \right) \right] \label{eq.db_dy.derive2} \end{eqnarray}\] \[\begin{eqnarray} \PartialDiff{c}{y_j} &=& \frac{1}{N\Delta} \left[ \left( <x^3> <x^1> - <x^2>^2 \right) x_j^2 \right. \nonumber \\ & & -\left( <x^4> <x^1> - <x^3> <x^2> \right) x_j \nonumber \\ & & \left. +\left( <x^4> <x^2> - <x^3>^2 \right) \right] \label{eq.dc_dy.derive2} \end{eqnarray}\] を得る。\(\sigma_{yj}\)の推定値としては \[\begin{equation} \sigma_{yj}=\sqrt{\frac{\sum_{i=1}^{N}e_i^2}{N-3}} \label{eq.sigma_yj} \end{equation}\] を使用する。ここで\(e_i\)はデータとフィッティング曲線との残差 \[\begin{equation} e_i = y_i-(ax_i^2+bx_i+c) = y_i-ax_i^2-bx_i-c \label{eq.ei} \end{equation}\] である。 (\ref{eq.da_dy.derive2})-(\ref{eq.ei})を (\ref{eq.sigma_a.original})-(\ref{eq.sigma_c.original})に代入して \[\begin{eqnarray} \sigma_a^2 &=& \sum_{j=1}^N \frac{1}{N^2\Delta^2} \left[ \left( <x^2> <x^0> - <x^1>^2 \right) x_j^2 \right. \nonumber \\ & & -\left( <x^3> <x^0> - <x^2> <x^1> \right) x_j \nonumber \\ & & \left. +\left( <x^3> <x^1> - <x^2>^2 \right) \right]^2 \nonumber \\ & & \frac{\sum_{i=1}^{N}e_i^2}{N-3} \nonumber \\ &=& \frac{1}{N^2(N-3)\Delta^2} \nonumber \\ & & \sum_{j=1}^N \left[ \left( <x^2> <x^0> - <x^1>^2 \right) x_j^2 \right. \nonumber \\ & & -\left( <x^3> <x^0> - <x^2> <x^1> \right) x_j \nonumber \\ & & \left. +\left( <x^3> <x^1> - <x^2>^2 \right) \right]^2 \nonumber \\ & & \sum_{i=1}^{N}\left(y_i-ax_i^2-bx_i-c\right)^2 \label{eq.sigma_a} \end{eqnarray}\] \[\begin{eqnarray} \sigma_b^2 &=& \sum_{j=1}^N \frac{1}{N^2\Delta^2} \left[ -\left( <x^3> <x^0> - <x^2> <x^1> \right) x_j^2 \right. \nonumber \\ & & +\left( <x^4> <x^0> - <x^2>^2 \right) x_j \nonumber \\ & & \left. -\left( <x^4> <x^1> - <x^3> <x^2> \right) \right]^2 \nonumber \\ & & \frac{\sum_{i=1}^{N}e_i^2}{N-3} \nonumber \\ &=& \frac{1}{N^2(N-3)\Delta^2} \nonumber \\ & & \sum_{j=1}^N \left[ -\left( <x^3> <x^0> - <x^2> <x^1> \right) x_j^2 \right. \nonumber \\ & & +\left( <x^4> <x^0> - <x^2>^2 \right) x_j \nonumber \\ & & \left. -\left( <x^4> <x^1> - <x^3> <x^2> \right) \right]^2 \nonumber \\ & & \sum_{i=1}^{N}\left(y_i-ax_i^2-bx_i-c\right)^2 \label{eq.sigma_b} \end{eqnarray}\] \[\begin{eqnarray} \sigma_c^2 &=& \sum_{j=1}^N \frac{1}{N^2\Delta^2} \left[ \left( <x^3> <x^1> - <x^2>^2 \right) x_j^2 \right. \nonumber \\ & & -\left( <x^4> <x^1> - <x^3> <x^2> \right) x_j \nonumber \\ & & \left. +\left( <x^4> <x^2> - <x^3>^2 \right) \right]^2 \nonumber \\ & & \frac{\sum_{i=1}^{N}e_i^2}{N-3} \nonumber \\ &=& \frac{1}{N^2(N-3)\Delta^2} \nonumber \\ & & \sum_{j=1}^N \left[ \left( <x^3> <x^1> - <x^2>^2 \right) x_j^2 \right. \nonumber \\ & & -\left( <x^4> <x^1> - <x^3> <x^2> \right) x_j \nonumber \\ & & \left. +\left( <x^4> <x^2> - <x^3>^2 \right) \right]^2 \nonumber \\ & & \sum_{i=1}^{N}\left(y_i-ax_i^2-bx_i-c\right)^2 \label{eq.sigma_c} \end{eqnarray}\] を得る。
No error is assumed for \(x_i\), and the error of \(y_i\) is represented by \(\sigma_{yi}\). Then the error propagation law tells that the estimation errors for \(a\), \(b\), and \(c\) are given by eqs (\ref{eq.sigma_a.original})-(\ref{eq.sigma_c.original}), respectively. From eqs (\ref{eq.a})-(\ref{eq.c}), (\ref{eq.da_dy.derive1})-(\ref{eq.dc_dy.derive1}) are obtained. Inserting eqs (\ref{eq.dx2y_dy})-(\ref{eq.dy_dy}) into these equations results in (\ref{eq.da_dy.derive2})-(\ref{eq.dc_dy.derive2}). For \(\sigma_{yi}\), eq. (\ref{eq.sigma_yj}) is assumed, where \(e_i\) is the residual between the data and fitting line (eq. \ref{eq.ei}). Inserting eqs. (\ref{eq.da_dy.derive2})-(\ref{eq.ei}) into (\ref{eq.sigma_a.original})-(\ref{eq.sigma_c.original}) results in (\ref{eq.sigma_a})-(\ref{eq.sigma_c}).
なお、\(N=3\)のときには誤差の推定は不可能であるので \(\sigma_a=\sigma_b=\sigma_c=0\)とする。
In case of \(N=3\), the error estimation is impossible; in this case, \(\sigma_a=\sigma_b=\sigma_c=0\) is returned.