計算に用いる時系列データを表す構造体。
A structure to represent the time series data
used for the computation.
◆戻り値(Return value)
傾き\(a\)、切片\(b\)、それらの推定誤差\(\sigma_a\), \(\sigma_b\)を並べた配列。
配列要素の順番は傾き、切片、傾きの推定誤差、切片の推定誤差の順とする。
切片は時系列データの定義域に関わらず直線の時刻\(t=0\)での値とする。
An array composed of the slope \(a\), intercept \(b\),
and their estimation errors \(\sigma_a\) and \(\sigma_b\).
The order of the array components is as follows:
the slope, the intercept, the estimation error of the slope,
and the estimation error of the intercept.
The intercept is defined as the value of the straight line at time \(t=0\)
regardless of the definition range of the time series data.
この例ではtrend[0]が傾き、trend[1]が切片、
trend[2]が傾きの推定誤差、trend[3]が切片の推定誤差となる。
In this example, trend[0] and trend[1] are
the slope and intercept,
and trend[2] and trend[3] are
the estimation erros of the slope and intercept, respectively.
◆計算式(Formula)
以下では最終的にプログラムで用いている式を赤色で示す。
Below, red colors show
the equations that are used in the program.
時系列データが時刻\(t=t_0,t_1,\cdots,t_{N-1}\)において定義されているものとし、
それらの時刻での値を\(f_0,f_1,\cdots,f_{N-1}\)とする。
この時系列データを最も良く説明する直線を\(f=at+b\)として
最小二乗フィッティングによって傾き\(a\)と切片\(b\)を求める。
観測方程式は
\[\begin{equation}
f_i=at_i+b \hspace{1em} (i=0,1,\cdots,N-1)
\label{eq.solve}
\end{equation}\]
であり、行列形式で
\[\begin{equation}
\begin{pmatrix}
f_0 \\
\vdots \\
f_{N-1}
\end{pmatrix}
=
\begin{pmatrix}
t_0 & 1 \\
\vdots & \vdots \\
t_{N-1} & 1
\end{pmatrix}
\begin{pmatrix}
a \\
b
\end{pmatrix}
\label{eq.solve.matrix}
\end{equation}\]
と書ける。この最小二乗解は
\[\begin{equation}
\begin{pmatrix}
t_0 & \cdots & t_{N-1} \\
1 & \cdots & 1
\end{pmatrix}
\begin{pmatrix}
f_0 \\
\vdots \\
f_{N-1}
\end{pmatrix}
=
\begin{pmatrix}
t_0 & \cdots & t_{N-1} \\
1 & \cdots & 1
\end{pmatrix}
\begin{pmatrix}
t_0 & 1 \\
\vdots & \vdots \\
t_{N-1} & 1
\end{pmatrix}
\begin{pmatrix}
a \\
b
\end{pmatrix}
\label{eq.solve.matrix.leastsquare}
\end{equation}\]
を解くことによって求められる。
左辺・右辺の係数行列を計算すると
\[\begin{equation}
\begin{pmatrix}
\sum_{i=0}^{N-1} t_i f_i \\
\sum_{i=0}^{N-1} f_i
\end{pmatrix}
=
\begin{pmatrix}
\sum_{i=0}^{N-1} t_i^2 & \sum_{i=0}^{N-1} t_i \\
\sum_{i=0}^{N-1} t_i & N
\end{pmatrix}
\begin{pmatrix}
a \\
b
\end{pmatrix}
\label{eq.solve.matrix.leastsquare.arrange}
\end{equation}\]
となり、クラメルの公式により
\[\begin{equation}
\textcolor[named]{Red}{
a=\frac{\Delta_a}{\Delta}
}
\label{eq.a}
\end{equation}\]
\[\begin{equation}
\textcolor[named]{Red}{
b=\frac{\Delta_b}{\Delta}
}
\label{eq.b}
\end{equation}\]
\[\begin{equation}
\textcolor[named]{Red}{\Delta}=
\begin{vmatrix}
\sum_{i=0}^{N-1} t_i^2 & \sum_{i=0}^{N-1} t_i \\
\sum_{i=0}^{N-1} t_i & N
\end{vmatrix}
\textcolor[named]{Red}{
=N\sum_{i=0}^{N-1} t_i^2-\left(\sum_{i=0}^{N-1} t_i\right)^2
}
\label{eq.Delta}
\end{equation}\]
\[\begin{equation}
\textcolor[named]{Red}{\Delta_a}=
\begin{vmatrix}
\sum_{i=0}^{N-1} t_i f_i & \sum_{i=0}^{N-1} t_i \\
\sum_{i=0}^{N-1} f_i & N
\end{vmatrix}
\textcolor[named]{Red}{
=N\sum_{i=0}^{N-1} t_i f_i
-\left(\sum_{i=0}^{N-1} t_i\right)\left(\sum_{i=0}^{N-1} f_i\right)
}
\label{eq.Delta_a}
\end{equation}\]
\[\begin{equation}
\textcolor[named]{Red}{\Delta_b}=
\begin{vmatrix}
\sum_{i=0}^{N-1} t_i^2 & \sum_{i=0}^{N-1} t_i f_i \\
\sum_{i=0}^{N-1} t_i & \sum_{i=0}^{N-1} f_i
\end{vmatrix}
\textcolor[named]{Red}{
=\left(\sum_{i=0}^{N-1} t_i^2\right)
\left(\sum_{i=0}^{N-1} f_i\right)
-\left(\sum_{i=0}^{N-1} t_i\right)
\left(\sum_{i=0}^{N-1} t_i f_i\right)
}
\label{eq.Delta_b}
\end{equation}\]
を得る。
Let the time series data be defined at time
\(t=t_0,t_1,\cdots,t_{N-1}\),
and the values of the data at these time samples be
\(f_0,f_1,\cdots,f_{N-1}\).
Let the straight line that best explains this time series data be \(f=at+b\).
We solve for the slope \(a\) and the intersect \(b\) of this straight line
by the least squares fitting.
The observation equation is eq. (\ref{eq.solve}),
which can be written in a matrix form as (\ref{eq.solve.matrix}).
The least squares solution of this equation can be obtained
by solving eq. (\ref{eq.solve.matrix.leastsquare}),
which can be arranged as (\ref{eq.solve.matrix.leastsquare.arrange}).
According to Cramer's formula, the solution of this equation is given by
eqs. (\ref{eq.a})-(\ref{eq.Delta_b}).
推定誤差の計算においては時刻\(t_i\)は誤差を持たないものとし、
各時刻での時系列データの値\(f_i\)は共通の誤差レベル
\[\begin{equation}
\sigma_f=\sqrt{\frac{\sum_{i=0}^{N-1}e_i^2}{N-2}}
\label{eq.sigma_f.definition}
\end{equation}\]
を持つものとする。ここで
\[\begin{equation}
e_i=f_i-(at_i+b)
\label{eq.residual}
\end{equation}\]
は時系列データとフィッティング直線の残差である。
(\ref{eq.residual})を(\ref{eq.sigma_f.definition})に代入して
\[\begin{eqnarray}
\textcolor[named]{Red}{\sigma_f^2}
&=& \frac{\sum_{i=0}^{N-1}(f_i-at_i-b)^2}{N-2}
\nonumber \\
&=& \frac{1}{N-2}
\sum_{i=0}^{N-1}(f_i^2+a^2t_i^2+b^2-2at_if_i-2bf_i+2abt_i)
\nonumber \\
&\textcolor[named]{Red}{=}&
\textcolor[named]{Red}{
\frac{1}{N-2}
\left(
\sum_{i=0}^{N-1}f_i^2
+a^2\sum_{i=0}^{N-1}t_i^2
+Nb^2
-2a\sum_{i=0}^{N-1}t_if_i
-2b\sum_{i=0}^{N-1}f_i
+2ab\sum_{i=0}^{N-1}t_i
\right)
}
\label{eq.sigma_f.result}
\end{eqnarray}\]
を得る。
なお\(N=2\)のときは誤差を推定しようが無いので\(\sigma_f=0\)とする。
誤差伝搬則より傾き\(a\)と切片\(b\)の推定誤差はそれぞれ
\[\begin{equation}
\sigma_a
=\sqrt{\sum_{j=0}^{N-1}\left(\PartialDiff{a}{f_j}\right)^2\sigma_f^2}
\label{eq.sigma_a.definition}
\end{equation}\]
\[\begin{equation}
\sigma_b
=\sqrt{\sum_{j=0}^{N-1}\left(\PartialDiff{b}{f_j}\right)^2\sigma_f^2}
\label{eq.sigma_b.definition}
\end{equation}\]
と書ける。(\ref{eq.a})(\ref{eq.b})式および
\(\Delta\)が\(f_j\)に依存しないこと(\ref{eq.Delta}式)より
\[\begin{equation}
\PartialDiff{a}{f_j}
=\PartialDiff{}{f_j}\frac{\Delta_a}{\Delta}
=\frac{1}{\Delta}\PartialDiff{\Delta_a}{f_j}
\label{eq.da_dfj}
\end{equation}\]
\[\begin{equation}
\PartialDiff{b}{f_j}
=\PartialDiff{}{f_j}\frac{\Delta_b}{\Delta}
=\frac{1}{\Delta}\PartialDiff{\Delta_b}{f_j}
\label{eq.db_dfj}
\end{equation}\]
であるので、
これらを(\ref{eq.sigma_a.definition})(\ref{eq.sigma_b.definition})式に代入して
\[\begin{equation}
\sigma_a
=\sqrt{\sum_{j=0}^{N-1}
\left(\frac{1}{\Delta}\PartialDiff{\Delta_a}{f_j}\right)^2\sigma_f^2}
=\frac{\sigma_f}{\Delta}
\sqrt{\sum_{j=0}^{N-1}\left(\PartialDiff{\Delta_a}{f_j}\right)^2}
\label{eq.sigma_a.use_Delta}
\end{equation}\]
\[\begin{equation}
\sigma_b
=\sqrt{\sum_{j=0}^{N-1}
\left(\frac{1}{\Delta}\PartialDiff{\Delta_b}{f_j}\right)^2\sigma_f^2}
=\frac{\sigma_f}{\Delta}
\sqrt{\sum_{j=0}^{N-1}\left(\PartialDiff{\Delta_b}{f_j}\right)^2}
\label{eq.sigma_b.use_Delta}
\end{equation}\]
を得る。ここで\(\Delta\)と\(\sigma_f\)が\(j\)によらない共通の値であることを用いた。
(\ref{eq.Delta_a})(\ref{eq.Delta_b})式より
\[\begin{eqnarray}
\PartialDiff{\Delta_a}{f_j}
&=& \PartialDiff{}{f_j}\left[
N\sum_{i=0}^{N-1} t_i f_i
-\left(\sum_{i=0}^{N-1} t_i\right)\left(\sum_{i=0}^{N-1} f_i\right)
\right]
\nonumber \\
&=& N\sum_{i=0}^{N-1} t_i\PartialDiff{f_i}{f_j}
-\left(\sum_{i=0}^{N-1} t_i\right)
\left(\sum_{i=0}^{N-1} \PartialDiff{f_i}{f_j}\right)
\nonumber \\
&=& N\sum_{i=0}^{N-1} t_i\delta_{ij}
-\left(\sum_{i=0}^{N-1} t_i\right)
\left(\sum_{i=0}^{N-1} \delta_{ij}\right)
\nonumber \\
&=& Nt_j-\sum_{i=0}^{N-1} t_i
\label{eq.Delta_a_diff}
\end{eqnarray}\]
\[\begin{eqnarray}
\PartialDiff{\Delta_b}{f_j}
&=& \PartialDiff{}{f_j}\left[
\left(\sum_{i=0}^{N-1} t_i^2\right)
\left(\sum_{i=0}^{N-1} f_i\right)
-\left(\sum_{i=0}^{N-1} t_i\right)
\left(\sum_{i=0}^{N-1} t_i f_i\right)
\right]
\nonumber \\
&=& \left(\sum_{i=0}^{N-1} t_i^2\right)
\left(\sum_{i=0}^{N-1} \PartialDiff{f_i}{f_j}\right)
-\left(\sum_{i=0}^{N-1} t_i\right)
\left(\sum_{i=0}^{N-1} t_i \PartialDiff{f_i}{f_j}\right)
\nonumber \\
&=& \left(\sum_{i=0}^{N-1} t_i^2\right)
\left(\sum_{i=0}^{N-1} \delta_{ij}\right)
-\left(\sum_{i=0}^{N-1} t_i\right)
\left(\sum_{i=0}^{N-1} t_i \delta_{ij}\right)
\nonumber \\
&=& \sum_{i=0}^{N-1} t_i^2-t_j\sum_{i=0}^{N-1} t_i
\label{eq.Delta_b_diff}
\end{eqnarray}\]
であるので
\[\begin{eqnarray}
\sum_{j=0}^{N-1}\left(\PartialDiff{\Delta_a}{f_j}\right)^2
&=& \sum_{j=0}^{N-1}\left[
Nt_j-\sum_{i=0}^{N-1} t_i
\right]^2
\nonumber \\
&=& \sum_{j=0}^{N-1}\left[
N^2 t_j^2
-2Nt_j\left(\sum_{i=0}^{N-1} t_i\right)
+\left(\sum_{i=0}^{N-1} t_i\right)^2
\right]
\nonumber \\
&=& N^2\sum_{j=0}^{N-1}t_j^2
-2N\left(\sum_{j=0}^{N-1}t_j\right)\left(\sum_{i=0}^{N-1} t_i\right)
+N\left(\sum_{i=0}^{N-1} t_i\right)^2
\nonumber \\
&=& N^2\sum_{i=0}^{N-1}t_i^2
-N\left(\sum_{i=0}^{N-1} t_i\right)^2
\nonumber \\
&=& N\left[N\sum_{i=0}^{N-1}t_i^2-\left(\sum_{i=0}^{N-1} t_i\right)^2\right]
\nonumber \\
&=& N\Delta
\label{eq.Delta_a_diff_sqsum}
\end{eqnarray}\]
\[\begin{eqnarray}
\sum_{j=0}^{N-1}\left(\PartialDiff{\Delta_b}{f_j}\right)^2
&=& \sum_{j=0}^{N-1}\left[
\sum_{i=0}^{N-1} t_i^2-t_j\sum_{i=0}^{N-1} t_i
\right]^2
\nonumber \\
&=& \sum_{j=0}^{N-1}\left[
\left(\sum_{i=0}^{N-1} t_i^2\right)^2
-2t_j\left(\sum_{i=0}^{N-1} t_i\right)\left(\sum_{i=0}^{N-1} t_i^2\right)
+t_j^2\left(\sum_{i=0}^{N-1} t_i\right)^2
\right]
\nonumber \\
&=& N\left(\sum_{i=0}^{N-1} t_i^2\right)^2
-2\left(\sum_{j=0}^{N-1}t_j\right)\left(\sum_{i=0}^{N-1} t_i\right)
\left(\sum_{i=0}^{N-1} t_i^2\right)
\nonumber \\
& & +\left(\sum_{j=0}^{N-1}t_j^2\right)\left(\sum_{i=0}^{N-1} t_i\right)^2
\nonumber \\
&=& N\left(\sum_{i=0}^{N-1} t_i^2\right)^2
-\left(\sum_{i=0}^{N-1} t_i\right)^2\left(\sum_{i=0}^{N-1} t_i^2\right)
\nonumber \\
&=& \left(\sum_{i=0}^{N-1} t_i^2\right)
\left[N\left(\sum_{i=0}^{N-1} t_i^2\right)
-\left(\sum_{i=0}^{N-1} t_i\right)^2\right]
\nonumber \\
&=& \left(\sum_{i=0}^{N-1} t_i^2\right)\Delta
\label{eq.Delta_b_diff_sqsum}
\end{eqnarray}\]
となり、これらを
(\ref{eq.sigma_a.use_Delta})(\ref{eq.sigma_b.use_Delta})
に代入して
\[\begin{equation}
\textcolor[named]{Red}{\sigma_a}
=\frac{\sigma_f}{\Delta}\sqrt{N\Delta}
\textcolor[named]{Red}{
=\sigma_f\sqrt{\frac{N}{\Delta}}
}
\label{eq.sigma_a.result}
\end{equation}\]
\[\begin{equation}
\textcolor[named]{Red}{\sigma_b}
=\frac{\sigma_f}{\Delta}\sqrt{\left(\sum_{i=0}^{N-1} t_i^2\right)\Delta}
\textcolor[named]{Red}{
=\sigma_f\sqrt{\frac{1}{\Delta}\sum_{i=0}^{N-1} t_i^2}
}
\label{eq.sigma_b.result}
\end{equation}\]
を得る。
For the computation of estimation errors,
times \(t_i\) are supposed to have no errors,
and the values \(f_i\) of the time series data at individual time samples
are supposed to have a common error level
given by Eq. (\ref{eq.sigma_f.definition}),
where \(e_i\) is the residual between the time series data
and the fitting line (Eq. \ref{eq.residual}).
Inserting Eq. (\ref{eq.residual}) into (\ref{eq.sigma_f.definition})
results in (\ref{eq.sigma_f.result});
in case of \(N=2\), the errors cannot be estimated,
and \(\sigma_f=0\) is used in this case.
According to the error propagation law,
estimation errors for the slope \(a\) and intercept \(b\)
are given by Eqs. (\ref{eq.sigma_a.definition})
and (\ref{eq.sigma_b.definition}), respectively.
From Eqs. (\ref{eq.a}) and (\ref{eq.b})
and the no dependency of \(\Delta\) on \(f_j\) (Eq. \ref{eq.Delta}),
we have Eqs. (\ref{eq.da_dfj}) and (\ref{eq.db_dfj}).
Inserting these equations into
(\ref{eq.sigma_a.definition}) and (\ref{eq.sigma_b.definition})
results in (\ref{eq.sigma_a.use_Delta})
and (\ref{eq.sigma_b.use_Delta}),
where the no dependencies of \(\Delta\) and \(\sigma_f\) on \(j\)
were used.
From Eqs. (\ref{eq.Delta_a}) and (\ref{eq.Delta_b}),
we have Eqs. (\ref{eq.Delta_a_diff}) and (\ref{eq.Delta_b_diff}),
and thus Eqs. (\ref{eq.Delta_a_diff_sqsum})
and (\ref{eq.Delta_b_diff_sqsum}).
Inserting these results into Eqs.
(\ref{eq.sigma_a.use_Delta}) and (\ref{eq.sigma_b.use_Delta}),
we finally obtain Eqs. (\ref{eq.sigma_a.result})
and (\ref{eq.sigma_b.result}).
(\ref{eq.Delta})-(\ref{eq.Delta_b}), (\ref{eq.sigma_f.result}),
(\ref{eq.sigma_b.result})式には\(t_i\)や\(t_i^2\)の和が登場し、
これらは時系列データのサンプリング間隔が一定であることを用いて解析式で計算できる。
時系列データのサンプリング間隔を\(\Delta t\)とおくと
\[\begin{equation}
t_i=t_0+i\Delta t
\label{eq.ti}
\end{equation}\]
であるので、
\[\begin{eqnarray}
\textcolor[named]{Red}{\sum_{i=0}^{N-1} t_i}
&=& \sum_{i=0}^{N-1}(t_0+i\Delta t)
\nonumber \\
&=& Nt_0+\Delta t\sum_{i=1}^{N-1}i
\nonumber \\
&\textcolor[named]{Red}{=}&
\textcolor[named]{Red}{
Nt_0+\Delta t\frac{N(N-1)}{2}
}
\label{eq.sumti}
\end{eqnarray}\]
\[\begin{eqnarray}
\textcolor[named]{Red}{\sum_{i=0}^{N-1} t_i^2}
&=& \sum_{i=0}^{N-1}(t_0+i\Delta t)^2
\nonumber \\
&=& \sum_{i=0}^{N-1}(t_0^2 + 2t_0 i \Delta t + i^2\Delta t^2)
\nonumber \\
&=& Nt_0^2 + 2t_0\Delta t\sum_{i=1}^{N-1}i + \Delta t^2\sum_{i=1}^{N-1}i^2
\nonumber \\
&\textcolor[named]{Red}{=}&
\textcolor[named]{Red}{
Nt_0^2 + t_0\Delta t N(N-1)
+\Delta t^2\frac{N(N-1)(2N-1)}{6}
}
\label{eq.sumti2}
\end{eqnarray}\]
である。
Eqs. (\ref{eq.Delta})-(\ref{eq.Delta_b}),
(\ref{eq.sigma_f.result}), and (\ref{eq.sigma_b.result})
have summations of \(t_i\) and \(t_i^2\),
which can be calculated analytically
because the time series data is sampled at a uniform interval.
The sample times are given by Eq. (\ref{eq.ti}),
where \(\Delta t\) is the sampling interval.
Therefore, we obtain Eqs. (\ref{eq.sumti}) and (\ref{eq.sumti2}).