高精度なインクリメンタル統計計算―分散値―

はじめに

別記事で，平均値の統計計算アルゴリズムの改善について軽く考察した．本記事ではその続きとして，分散の計算アルゴリズムを改良する．

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)

パラメータと関数の設定

実は平均のときにも，執筆途中で記事に載せているプログラムから仕様を変更したのだが，ここではパラメータと分散を計算する関数を外部手続きとして別ファイルに記述して，各プログラムで設定を使い回せるようにしている．

具体的には，次のようにパラメータ用のファイルと関数用のファイルを作成して，

stats_parameters.f90

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

stats_function.f90

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を使ってモジュールを読み込むようにして，

main_program.f90

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コードおよび実行結果は，

stats_variance_naive.f90

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で実装すると，

stats_variance_online.f90

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}

stats_variance_online_merge.f90

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: カハンの加算アルゴリズム) を加えてみる．

stats_variance_online_merge_kahan.f90

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)$を次のように定義する．

stats_function.f90

        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を加えない場合とほとんど結果に差異が生じないことがわかった．