はじめに
別記事で，平均値の統計計算アルゴリズムの改善について軽く考察した．本記事ではその続きとして，分散の計算アルゴリズムを改良する．
to-kei.netより，
統計学において、分散とは数値データのばらつき具合を表すための指標です。
また，分散$s^2$は，$n$を観測値の数，$x_1,x_2,\cdots,x_n$を一つ一つの観測値，$\bar{x}$をこれらの観測値の平均とすると，
$$
s^2 = \frac{1}{n}\sum_{i=1}^n(x_i - \bar{x})^2
$$
$i = 1,2,\cdots,n$を経過時間だと考えると，上記は離散時間$i$に対して分散を計算しているともいえる．
これ以後，アルゴリズム確認はFortranで，そのコンパイラは
$ gfortran -v
...
gcc version 7.4.0 (Ubuntu 7.4.0-1ubuntu1~18.04.1)
パラメータと関数の設定
実は平均のときにも，執筆途中で記事に載せているプログラムから仕様を変更したのだが，ここではパラメータと分散を計算する関数を外部手続きとして別ファイルに記述して，各プログラムで設定を使い回せるようにしている．
具体的には，次のようにパラメータ用のファイルと関数用のファイルを作成して，
module stats_parameters
use,intrinsic :: iso_fortran_env
implicit none
integer(int8 ),parameter :: kp = selected_real_kind(6) !6=single precision, 15=double precision
integer(int32),parameter :: N = 10000000 !データ数
integer(int32),parameter :: DT = 1 !離散時間間隔
integer(int32),parameter :: M_DT = 50000 !マージ時間間隔
end module stats_parameters
module stats_function
use,intrinsic :: iso_fortran_env
implicit none
contains
function f(t) result(y)
use,intrinsic :: iso_fortran_env
use stats_parameters
implicit none
real(kp) x,y
integer(int32) t
y = 0.01_kp !サンプル関数1: 定数関数
! x = 2*acos(-1.0_kp)/(N/DT)*t
! y = 1.0_kp+sin(x) !サンプル関数2: 三角関数
end function f
end module stats_function
メインプログラムのほうでuseを使ってモジュールを読み込むようにして，
program main
use,intrinsic :: iso_fortran_env
use stats_parameters
use stats_function
implicit none
...
コンパイル時に関連するファイルをまとめてコンパイルする (順番に注意)．
$ gfortran stats_parameters.f90 stats_function.f90 main_program.f90
不偏標本分散
ここで計算している分散は不偏標本分散 (unbiased variance) である．超簡単に説明すると，標本分散 (sample variance) とは上記はじめにで示した定義式の$1/n$をそのまま用いて，合計$n$個のデータ (これを標本という) における問題限定の値として分散を扱う．一方で不偏分散は，$1/n$を$1/(n-1)$とすることで，標本を含めた全体の枠組み (母集団) における問題の値として分散を扱う．例えば，学校で学力テストを行って，学校単位で分散 (から出す偏差値) を考える場合は標本分散，全国の学校全体で分散を考える場合は不偏分散を利用するようなものである．詳細を知りたければ，to-kei.net等を確認する．
単純な例題 (離散時間間隔DT = 1)
平均のときと同様，もっとも単純な問題の一つとして，ここでは時間関数を
$$
f(t) = 0.01,\quad 1\le t\le N
$$
とする．すなわち，単精度実数に対して，(時間的には不変となっているが) 時間関数として常に0.010を取る関数を考える．この関数の平均値は解析的には任意のNに対して0.01であり，分散値は0である．N = 10,000,000，離散時間間隔DT = 1で離散化した関数$f$に対してまずは考えていく．
Naïve algorithm
Qiita記事：分散 = 2乗の平均 - 平均の2乗やWikipedia: Naïve algorithmで紹介されている，Naïve algorithmをまずは実装してみる．(本記事ではインクリメンタルなアルゴリズムの改良を目指すため，最初の定義に従った「平均を求めてから分散を求める」実直な計算は省略する．) そのFortranコードおよび実行結果は，
program stats_variance_naive
use,intrinsic :: iso_fortran_env
use stats_parameters
use stats_function
implicit none
real(kp) sum,sum_sq
integer(int32) t
sum = 0.0_kp
sum_sq = 0.0_kp
print*,'Naïve algorithm'
do t=1,N,DT
sum = sum+f(t)
sum_sq = sum_sq+f(t)*f(t)
enddo
print*,(sum_sq-(sum*sum)/(N/DT))/(N/DT-1) !unbiased variance
end program stats_variance_naive
$ gfortran stats_parameters.f90 stats_function.f90 stats_variance_naive.f90
$ ./a.out
Naïve algorithm
1.72239888E-05
平均のとき同様，表題にも書かれているDT (離散時間間隔) をわざわざパラメータとして与えているのは，そのうち役に立つので実装にだけ気をつけてひとまず脇においておく．この結果1.72239888E-05を改善していく (0に近づける)．
Online algorithm (Welford's online algorithm)
WikipediaではWelford's online algorithm，Qiita記事ではさらに誤差を抑えるアルゴリズムとして紹介されている漸近アルゴリズムを実装する．なお，Wikipediaにあるように，「sum of squares of differences from the current mean」として$M_{2,n}$を計算し，最後にデータ数 (不偏分散の場合はデータ数-1) で割って分散を出すという実装のほうが高精度である．Qiita記事では平均値についての式展開が載せられているため，こちらでは分散値を計算するための$M_{2,n}$の漸化式展開を挙げておく．これも平均値同様，合計を蓄積しないようにして丸め誤差を小さくする工夫のされた式になる．各文字定義は以下の表に従う (Wikipediaより)．
| Characters | Explanation |
|---|---|
| $x_i$ | 時刻$i$における値 |
| $\bar{x}n = (\sum{i=1}^n x_i)/n$ | $i=1-n$の平均値 |
| $M_{2,n} = \sum_{i=1}^n (x_i-\bar{x}_n)^2$ | $i=1-n$までの平均からの差の2乗和 |
\begin{align}
M_{2,n} &= \sum_{i=1}^n (x_i-\bar{x}_n)^2 \\
&= \sum_{i=1}^{n-1} (x_i-\bar{x}_n)^2 + (x_n-\bar{x}_n)^2 \\
&= \sum_{i=1}^{n-1} \left\{x_i-\left[\bar{x}_{n-1}+\frac{1}{n}(x_n-\bar{x}_{n-1})\right]\right\}^2 + (x_n-\bar{x}_n)^2\qquad\text{∵平均に対する漸近アルゴリズム} \\
&= \sum_{i=1}^{n-1} \left\{(x_i-\bar{x}_{n-1})-\frac{1}{n}(x_n-\bar{x}_{n-1})\right\}^2 + (x_n-\bar{x}_n)^2 \\
&= \sum_{i=1}^{n-1} (x_i-\bar{x}_{n-1})^2-\frac{2}{n}(x_n-\bar{x}_{n-1})\sum_{i=1}^{n-1} (x_i-\bar{x}_{n-1})+\frac{1}{n^2}(x_i-\bar{x}_{n-1})^2(n-1) + (x_n-\bar{x}_n)^2 \\
\end{align}
ここで，右辺第1項は$M_{2,n-1}$，第2項は$0$，第3・4項は
\begin{align}
& \frac{1}{n^2}(x_i-\bar{x}_{n-1})^2(n-1) + (x_n-\bar{x}_n)^2 \\
&\quad= \frac{1}{n^2}(x_i-\bar{x}_{n-1})(n-1)\left(\frac{n}{n-1}x_n-\frac{n}{n-1}\bar{x}_n\right)+(x_n-\bar{x}_n)\left(x_n-\frac{n-1}{n}\bar{x}_{n-1}-\frac{1}{n}x_n\right) \\
&\quad= \frac{1}{n}(x_n-\bar{x}_{n-1})(x_n-\bar{x}_n)+(x_n-\bar{x}_n)\frac{n-1}{n}(x_n-\bar{x}_{n-1}) \\
&\quad= (x_n-\bar{x}_n)(x_n-\bar{x}_{n-1})
\end{align}
以上より，
M_{2,n} = M_{2,n-1}+(x_n-\bar{x}_n)(x_n-\bar{x}_{n-1})
これをFortranで実装すると，
program stats_variance_online
use,intrinsic :: iso_fortran_env
use stats_parameters
use stats_function
implicit none
real(kp) avg,avg_new,sum_online
integer(int32) t
avg = 0.0_kp
avg_new = 0.0_kp
sum_online = 0.0_kp
print*,'Welford''s online algorithm' !use M2,n
do t=1,N
if(mod(t,DT)==0)then
avg_new = avg*(1.0_kp-real(DT,kp)/t)+f(t)*real(DT,kp)/t !online for average (mean)
sum_online = sum_online+(f(t)-avg)*(f(t)-avg_new)
avg = avg_new
endif
enddo
print*,sum_online/(N/DT-1) !unbiased variance
end program stats_variance_online
$ gfortran stats_parameters.f90 stats_function.f90 stats_variance_online.f90
$ ./a.out
Welford's online algorithm
2.14895586E-08
先程の1.72239888E-05に比べて3桁だけオーダーが小さくなってくれた．
Online+Parallel algorithm (Online+Merge algorithm)
次に，このコードにMergeのアルゴリズム (Qiita記事ではマージする，WikipediaではParallel algorithm) を追加する．平均の場合はQiita記事に式展開がある (Wikipediaにはアルゴリズムに使える式のみ)．文字定義は先の表に加えて，以下の表に倣う (Wikipedia準拠)．
| Characters | Explanation |
|---|---|
| $A,B$ | 対象とする全データを$A,B$の二つに分割 |
| $\delta$ | $\bar{x}_B-\bar{x}_A$ |
| $n_X$ | 対象とする全データ数: $n_A+n_B$ |
\begin{align}
M_{2,X} &= \sum_{i=1}^{n_X}(x_i-\bar{x}_{X})^2 \\
&= \sum_{i=1}^{n_A}(x_i-\bar{x}_{X})^2+\sum_{i=n_A+1}^{n_X}(x_i-\bar{x}_{X})^2 \\
&= \sum_{i=1}^{n_A}\left[x_i-\left(\bar{x}_{A}+\delta\frac{n_B}{n_X}\right)\right]^2+\sum_{i=n_A+1}^{n_X}\left[x_i-\left(\bar{x}_{B}-\delta\frac{n_A}{n_X}\right)\right]^2\qquad\text{(∵$\bar{x}_X = (n_A\bar{x}_A+n_B\bar{x}_B)/n_X$)} \\
&= \sum_{i=1}^{n_A}\left[(x_i-\bar{x}_{A})^2-2(x_i-\bar{x}_{A})\delta\frac{n_B}{n_X}+\delta^2\frac{n_B^2}{n_X^2}\right]+\sum_{i=n_A+1}^{n_X}\left[(x_i-\bar{x}_{B})^2+2(x_i-\bar{x}_{B})\delta\frac{n_A}{n_X}+\delta^2\frac{n_A^2}{n_X^2}\right] \\
&= M_{2,A}-2\delta\frac{n_B}{n_X}\sum_{i=1}^{n_A}(x_i-\bar{x}_A)+\delta^2\frac{n_An_B^2}{n_X^2}+M_{2,B}+2\delta\frac{n_A}{n_X}\sum_{i=n_A+1}^{n_X}(x_i-\bar{x}_B)+\delta^2\frac{n_A^2n_B}{n_X^2} \\
&= M_{2,A}+M_{2,B}-2\delta\frac{1}{n_X}\left(n_B\sum_{i=1}^{n_A}x_i-n_A\sum_{i=n_A+1}^{n_X}x_i\right)+2\delta\frac{n_An_B}{n_X}\bar{x}_A-2\delta\frac{n_An_B}{n_X}\bar{x}_B+\delta^2\frac{n_An_B(n_A+n_B)}{n_X^2} \\
&= M_{2,A}+M_{2,B}-2\delta\frac{1}{n_X}(n_An_B\bar{x}_A-n_An_B\bar{x}_B)-\delta^2\frac{n_An_B}{n_X}\qquad\text{(∵$n_X = n_A+n_B$)} \\
&= M_{2,A}+M_{2,B}+2\delta^2\frac{n_An_B}{n_X}-\delta^2\frac{n_An_B}{n_X} \\
&= M_{2,A}+M_{2,B}+\delta^2\frac{n_An_B}{n_X}
\end{align}
program stats_variance_online_merge
use,intrinsic :: iso_fortran_env
use stats_parameters
use stats_function
implicit none
real(kp) avg_part,avg_new,avg_merge
real(kp) sum_part,sum_merge
integer(int32) t,t_part,ct
avg_part = 0.0_kp
avg_merge = 0.0_kp
sum_part = 0.0_kp
sum_merge = 0.0_kp
ct = 0
print*,'Online+Parallel algorithm (Online+Merge algorithm)'
do t=1,N
if(mod(t,DT)==0)then
t_part = t-ct*M_DT !ここでマージごとに次式の分母をリセットできるので和算のオーダーが離れづらくなる
avg_new = avg_part*(1.0_kp-real(DT,kp)/t_part)+f(t)*real(DT,kp)/t_part !online for average (mean)
sum_part = sum_part+(f(t)-avg_new)*(f(t)-avg_part) !online for variance
avg_part = avg_new
endif
if(mod(t,M_DT)==0)then
ct = t/M_DT
!パートAが更新前の_merge，データ量はt-M_DT．パートBが_part，データ量はM_DT
sum_merge = sum_merge+sum_part+(avg_part-avg_merge)**2*(t-M_DT)*M_DT/t !merge for variance
avg_merge = 1.0_kp/ct*((ct-1)*avg_merge+avg_part) !merge for average
avg_part = 0.0_kp !マージごとに個別の平均値はリセット
sum_part = 0.0_kp !マージごとに個別のsumはリセット
endif
enddo
print*,sum_merge/(N/DT-1) !unbiased variance
end program stats_variance_online_merge
t = 50,000ごとに (M_DT = 50000) マージするとして，
$ gfortran stats_parameters.f90 stats_function.f90 stats_variance_online_merge.f90
$ ./a.out
Online+Parallel algorithm (Online+Merge algorithm)
5.86959051E-16
前回のOnlineアルゴリズム2.14895586E-08と比べても8桁ほど相対誤差を小さくすることができた．この理由は，平均の記事でMergeアルゴリズムを平均計算に対して実装した際，Onlineアルゴリズムのみの結果と比較して-1.7902774972853197E-002 → 9.8944689321811552E-006と4桁だけ小さくなったことから，今回はこのオーダーの2乗でアルゴリズムが効いたと考えられる．
Online+Parallel (Merge)+Kahan algorithm
さらに，Kahanのアルゴリズム (Wikipedia: カハンの加算アルゴリズム) を加えてみる．
program stats_variance_online_merge_kahan
use,intrinsic :: iso_fortran_env
use stats_parameters
use stats_function
implicit none
real(kp) avg_part,avg_new,avg_merge
real(kp) sum_part,sum_merge
real(kp) yk,tk,ck
integer(int32) t,t_part,ct
avg_part = 0.0_kp
avg_merge = 0.0_kp
sum_part = 0.0_kp
sum_merge = 0.0_kp
ck = 0.0_kp
ct = 0
print*,'Online+Parallel (Merge) +Kahan algorithm'
do t=1,N
if(mod(t,DT)==0)then
t_part = t-ct*M_DT !ここでマージごとに次式の分母をリセットできるので和算のオーダーが離れづらくなる
avg_new = avg_part*(1.0_kp-real(DT,kp)/t_part)+f(t)*real(DT,kp)/t_part !online for average (mean)
!sum = (sum-ck)+f(t)
! = sum+(f(t)-ck)
! = sum+yk
! = tk
yk = (f(t)-avg_new)*(f(t)-avg_part)-ck
tk = sum_part+yk
ck = (tk-sum_part)-yk
sum_part = tk
avg_part = avg_new
endif
if(mod(t,M_DT)==0)then
ct = t/M_DT
!パートAが更新前の_merge，データ量はt-M_DT．パートBが_part，データ量はM_DT
sum_merge = sum_merge+sum_part+(avg_part-avg_merge)**2*(t-M_DT)*M_DT/t !merge for variance
avg_merge = 1.0_kp/ct*((ct-1)*avg_merge+avg_part) !merge for average
avg_part = 0.0_kp !マージごとに個別の平均値はリセット
sum_part = 0.0_kp !マージごとに個別のsumはリセット
endif
enddo
print*,sum_merge/(N/DT-1) !unbiased variance
end program stats_variance_online_merge_kahan
$ gfortran stats_parameters.f90 stats_function.f90 stats_variance_online_merge_kahan.f90
$ ./a.out
Online+Parallel (Merge) +Kahan algorithm
5.87027555E-16
Kahan導入前の5.86959051E-16に比べて，ほんっの少しだけ悪くなった．予想はしていたが，平均のとき同様，Kahanは他のアルゴリズムとの相性が悪いっぽい．
一旦まとめ
以上をまとめると，
Naïve algorithm
1.72239888E-05
Welford's online algorithm
2.14895586E-08
Online+Parallel algorithm (Online+Merge algorithm)
5.86959051E-16
Online+Parallel (Merge) +Kahan algorithm
5.87027555E-16
やっぱりMergeが強くてKahanが微妙い．
単純な例題 (離散時間間隔DT = 1000)
平均のとき同様，離散時間間隔をDT = 1000にしてみる．計算精度は単精度のまま (パラメータファイルのkp = selected_real_kind(6))．
Naïve algorithm
-5.60339719E-10
Welford's online algorithm
6.42166051E-17
Online+Parallel algorithm (Online+Merge algorithm)
1.53357020E-15
Online+Parallel (Merge) +Kahan algorithm
1.53357020E-15
平均のときと違う．平均のときはMergeアルゴリズムを加えるとDT = 1000でも結果が改善されていたが，今回はOnlineアルゴリズム (漸近アルゴリズム) のみの方がゼロに近いので解析値に近い．理由として考えられるのは，平均のときには1/tが1/t_partとなって，マージするたびにt_partを初期化することにより結果の改善が図れたのが，分散のOnlineアルゴリズムには分母に1/tが存在しないため，特に改善の方向にマージが作用しない (DT = 1のときに改善されているのは平均の改善作用)．さらに，離散時間間隔が大きくなることで計算に使われる総データ数が少なくなり，Onlineアルゴリズムのみの場合は1/tが大きくなりすぎないことで結果が改善される (2.14895586E-08→6.42166051E-17)．しかし総データ数が少なくなると平均に対するMergeアルゴリズムの改善幅が小さくなり，結果としてDT = 1000の場合はOnlineのみに対してOnline+Mergeで改悪されたと考えられる．
倍精度演算
計算精度を倍精度に (パラメータファイルのkp = selected_real_kind(15)) にしてみる．DT = 1の場合，
Naïve algorithm
9.4421703017156724E-015
Welford's online algorithm
4.8115845668530339E-030
Online+Parallel algorithm (Online+Merge algorithm)
1.1024187308774392E-032
Online+Parallel (Merge) +Kahan algorithm
1.1024187308774392E-032
DT = 1000の場合，
Naïve algorithm
-3.7895701400596654E-017
Welford's online algorithm
6.3421798570560140E-034
Online+Parallel algorithm (Online+Merge algorithm)
3.2440905470277344E-032
Online+Parallel (Merge) +Kahan algorithm
3.2440905470277344E-032
まぁそんなもんでしょうねという感じ．特に傾向は単精度のときと変わらない．
CPU時間計測
参考程度に，ある基準時点からのCPU時間を取得する組み込みサブルーチンcpu_timeを使ってループの計算時間を計ってみる．計算条件は，倍精度 (kp = selected_real_kind(15))，N = 100,000,000，DT = 1，M_DT = 50000として，3回計算した平均をとると (単位は全て秒)，
| Algorithm | Naïve | Online | Online+Merge | Online+Merge+Kahan |
|---|---|---|---|---|
| Avg. Time [s] | 1.1875 | 2.614583333 | 2.609375 | 2.59375 |
まとめ
平均アルゴリズムに続いて，インクリメンタルな分散アルゴリズムの改良を試みた．結果は，Naïve algorithmに比べて，Welford's online algorithm (漸近アルゴリズム) にすることで大幅に改善された．Online+Parallel (Merge) algorithmは，総データ数 (離散時間間隔) に依存してOnlineのみの結果に対する改善・改悪が変化した．Kahan algorithmを加えると，少なくとも今回の検証では特に改善も改悪もされなかった．平均アルゴリズムとの一貫性を考えると，分散のインクリメンタルな計算には，Online+Parallel (Merge) algorithmが無難だと考えられる．
追加: 関数を三角関数にして検証
せっかく関数を変更しやすくしたので，三角関数でアルゴリズムを検証してみた．パラメータと関数の設定で示したように，関数f(t)を総データ数N/DTに対して1周期分$0-2\pi$となる$1+\sin(x)$を次のように定義する．
x = 2*acos(-1.0_kp)/(N/DT)*t
y = 1.0_kp+sin(x) !サンプル関数2: 三角関数
この関数の解析的な平均値は$1$，分散値は$0.5$である．データ数N = 10,000,000，マージ時間間隔M_DT = 50000は固定．ついでで平均値も出している．全て相対誤差で表示．
| Precision | DT | Relative error | Naïve | Online | Online+Merge | Online+Merge+Kahan |
|---|---|---|---|---|---|---|
| Single | 1 | Average | 1.0674E-02 | -9.3600E-02 | 2.4199E-05 | 2.4199E-05 |
| Variance | -1.4937E-02 | -1.4452E-01 | 6.4969E-05 | 6.4969E-05 | ||
| Single | 1000 | Average | -2.3842E-07 | -2.1696E-05 | 2.3842E-07 | 2.3842E-07 |
| Variance | -2.8491E-05 | -1.0109E-04 | -1.0014E-04 | -1.0014E-04 | ||
| Double | 1 | Average | 2.6557E-13 | 8.0824E-14 | -2.2204E-16 | -2.2204E-16 |
| Variance | -1.0000E-07 | -1.0000E-07 | -1.0000E-07 | -1.0000E-07 | ||
| Double | 1000 | Average | 4.4409E-16 | -1.7319E-14 | 2.2204E-16 | 2.2204E-16 |
| Variance | -1.0000E-04 | -1.0000E-04 | -1.0000E-04 | -1.0000E-04 |
考察
色々と言えそうだが特に注目したいのは，
- OnlineのみだとNaïveに対するアドバンテージがほぼない，むしろ悪くなっていることの方が多い
- Online+MergeだとNaïveに対して，データ数に依存する傾向はあるがほとんどの場合，結果は改善される
- Online+MergeとOnline+Merge+Kahanはまったく同じ結果に
あたりである．ここから得られる結論としては，Online+Mergeがもっとも精度の高いアルゴリズムなので採用すべきであり，これは単純な問題として関数を$0.01$で与えた場合の結論に一致する．
三角関数での検証を踏まえたまとめ
まとめにおける，単純な問題として関数を$0.01$で与えた場合の検証に加え，三角関数でも同様にアルゴリズムの検証をした．もっとも精度の高いアルゴリズム (の組み合わせ) はOnline+Merge algorithmであり，Online algorithmのみの場合はNaïve algorithmより結果が悪くなる可能性が高く，Online+Merge+Kahan algorithmの場合はKahan algorithmを加えない場合とほとんど結果に差異が生じないことがわかった．