数値的に安定かつ高速な移動平均と移動分散の実装

本記事では、系列の各時点における直近 \(w\) 点の平均と分散を、次の 3 つの条件を同時に満たすように計算する方法を紹介します。（1）計算量が窓幅 \(w\) に依存せず \(O(T)\) であること、（2）累積和と要素ごとの演算だけで構成され、ベクトル化できること、（3）系列の水準やドリフトに左右されず、数値的に安定していることです。まず、移動和と移動平均を「チャンクごとの累積和」で安定して求めます（§1）。次に、同じチャンク構造を拡張して移動分散を求めます（§2）。用いるのは、Welford のオンラインアルゴリズム[1]と、それを一般化した Chan らのマージ式[2]という古典的な手法です。本記事は、これらを移動窓の計算に組み合わせる方法を整理したものです。PyTorch 向けの実装は、開発しているライブラリ qfeval-functions で msum・mvar として公開しています。

基本的な考え方は次の 3 点です。

移動和：系列を窓幅 \(w\) と同じ長さのチャンクに区切ると、どの窓も「前のチャンクの接尾辞＋次のチャンクの接頭辞」という最大 2 つの部分に分解できます。すべてのチャンクの接頭辞と接尾辞の和は、累積和を使って一括で計算できます。これらを組み合わせれば、移動和と移動平均を求められます（§1）。
移動分散への拡張（部分ごとの統計量）：各部分について保持する量を、「和」から「要素数・平均・偏差二乗和の 3 つ組」に置き換えます。このとき統計量は、その部分自身に含まれる要素（チャンクの先頭または末尾の値）を基準にした偏差から計算します。分散は値全体を平行移動しても変わらないため、計算結果を変えずに、誤差を窓内の局所的な値のスケールに抑えられます（§2.3）。
移動分散への拡張（統計量の結合）：2 つの部分の分散統計量を Chan らのマージ式で結合し、1 つの窓の統計量を求めます（§2.4）。

中核となるのは次のマージ式です。重なりのない 2 つの部分集合 \(A, B\) があり、それぞれの要素数 \(n\)・平均 \(\mu\)・偏差二乗和 \(M = \sum (x - \mu)^2\) が得られているとき、

\[ \boxed{\;M_{A\cup B} = M_A + M_B + \frac{n_A\,n_B}{n_A + n_B}\,\delta^2,\qquad \delta = \mu_B - \mu_A\;} \]

が成り立ちます（導出は §2.4）。移動分散は \(\mathrm{Var} = M_{A\cup B}/(w - \mathrm{ddof})\) です。この分解と結合をテンソル演算で表すと、すべての時点の移動分散を、4 本の累積和と要素ごとの演算だけで求められます。

\(\mathrm{ddof}\)（delta degrees of freedom）は、分散の分母を要素数よりいくつ小さくするかを指定する値です。記号は numpy や pandas の引数名に合わせています。標本分散（不偏分散）なら \(\mathrm{ddof} = 1\)、母分散なら \(0\) を使います（§2.4）。

対象とする問題

長さ \(T\) の系列 \(x_1, \dots, x_T\) と窓幅 \(w\) に対し、各時点 \(t\) について直近の窓 \((x_{t-w+1}, \dots, x_t)\) の和・平均・分散を求めます。出力は入力と同じ長さで、窓が系列からはみ出す先頭の \(w-1\) 時点は NaN とし、窓に NaN が含まれる時点も NaN とします（pandas の rolling(w) と同じ規約です）。

金融時系列には、単純な実装では精度や速度の問題が生じやすい条件が揃っています。例えば、価格の水準（\(10^4\)〜\(10^6\) など）は変動幅（\(1\) 程度）より桁違いに大きい一方で、GPU では有効桁が約 7 桁の float32 がよく使われます。さらに、数千系列、数十万時点のデータを一括で処理する場合もあります。そのため、次の 4 つの性質が必要です。

計算量が \(O(T)\)：計算時間が窓幅 \(w\) に依存しないこと。日次データでは、\(w = 250\) のように長い窓も一般的に使われます。
ベクトル化できる：全体がテンソル演算（累積和と要素ごとの演算）で構成され、時点方向に明示的な逐次ループを持たないこと。これにより、GPU や SIMD で効率よく実行できます。
数値的に安定している：誤差が窓内の局所的な変動のスケールで決まり、系列全体の水準・ドリフト・長さに依存しないこと。
NaN・±inf が正確に伝播する：欠損値や無限大が、それを含む窓だけに影響し、含まない窓の結果には影響しないこと。

最初の 2 つだけであれば単純な公式で、3 つ目だけであれば窓ごとに再計算する方法で満たせます。しかし、4 つすべてを同時に満たすには工夫が必要です。本記事では、まず単純な方法による移動和が長い系列で精度を失う理由を確認し、チャンク分割によって安定化します（§1）。その後、この構造を拡張して移動分散を構成します（§2）。

1. 移動和・移動平均：チャンク累積和

移動分散を扱う前に、まず移動和を考えます。移動平均は移動和を \(w\) で割るだけで求められるため、以下では移動和に絞って説明します。

1.1. 大域的な累積和の差分は系列長とともに精度が劣化する

移動和を \(O(T)\) で求める一般的な方法は、大域的な累積和の差分 \(S_t - S_{t-w}\)（\(S_t = \sum_{i \le t} x_i\)）を使うものです。しかし、この方法では系列が長くなるほど精度が低下します。水準 \(\mu\) の系列では \(S_t \approx t\mu\) と大きくなり、これを浮動小数点数として表現するだけで、絶対誤差が \(\varepsilon\, t\,|\mu|\) 程度になるためです。ここで、\(\varepsilon\) は丸め単位であり、float32 では約 \(6 \times 10^{-8}\) です。窓和の真値は \(w\mu\) 程度なので、相対誤差は \(\varepsilon\, t/w\) のオーダーとなり、系列の後半ほど大きくなります。float32 で \(t = 10^7\)、\(w = 100\) の場合、この誤差は 0.6% に達します。

この誤差は、累積和の計算アルゴリズムを工夫しても避けられません。逐次加算でもペアワイズ加算でも、計算結果の \(S_t\) を 1 つの浮動小数点数に丸めた時点で、その値には最良でも \(\varepsilon\,|S_t|\) 程度の表現誤差があります。\(S_t\) と \(S_{t-w}\) が同符号で値も近ければ、差分の演算自体は正確に行えます（Sterbenz の補題）。しかし、両者はすでに真値から \(\varepsilon\,t\,|\mu|\) 程度ずつずれているため、小さな真値 \(w\mu + O(w\sigma)\) に対して、その誤差が残ります。次節のチャンク分割は、累積和を最大 \(w\) 項で区切り、この誤差の蓄積を抑える方法とも解釈できます。

1.2. チャンク分割：どの窓も 2 つ以下の部分に分けられる

系列の先頭に NaN を \(w\) 個以上追加し、全体の長さを \(w\) の倍数に揃えます（個数の理由は §6 を参照）。そのうえで、先頭から順に長さ \(w\) のチャンクに区切ります。窓とチャンクの長さがどちらも \(w\) なので、どの窓も連続する最大 2 つのチャンクにまたがります。すなわち、チャンク \(c\) の中で終わる窓は、チャンク \(c-1\) の接尾辞（suffix）とチャンク \(c\) の接頭辞（prefix）を連結したものとして、

\[ \text{window} = \underbrace{\text{suffix of chunk } c-1}_{A\ (\text{length } w - k)}\ +\ \underbrace{\text{prefix of chunk } c}_{B\ (\text{length } k)} \qquad (1 \le k \le w) \]

と分解できます。\(k = w\) のとき、接頭辞はチャンク全体となり、接尾辞は空になります。各チャンクについて、長さ 1 から \(w\) までのすべての接頭辞の和は、先頭からの累積和で一括計算できます。同様に、すべての接尾辞の和は逆向きの累積和で一括計算できます。チャンクは \(T/w\) 個あり、1 本の累積和に必要な計算量はチャンクあたり \(O(w)\) なので、全体の計算量は \(O(T)\) です。

チャンク分割による窓の分解（w = 4）。チャンク 2 の中で終わる 4 つの窓が、いずれもチャンク 1 の接尾辞（青）とチャンク 2 の接頭辞（オレンジ色）に分かれ、チャンク境界で 1 回加算すると窓和になる。接尾辞の和は逆向きの累積和、接頭辞の和は累積和で一括計算できる。 — 図 1：チャンク分割による窓の分解（\(w = 4\)）どの窓も「前チャンクの接尾辞＋次チャンクの接頭辞」に分解できる。窓和は、接尾辞の和と接頭辞の和の 1 回の加算で求まる。各部分の和は、全チャンクで取った両向きの累積和（青とオレンジ色の矢印）が一括で与える。

1.3. 部分和から移動和・移動平均を求める

窓和は、接尾辞の和と接頭辞の和を 1 回加算するだけで求まります。各部分の和は最大 \(w\) 項からなるため、誤差は局所的な値のスケールに \(O(\varepsilon w)\) を掛けた範囲に収まります。大域累積和（§1.1）とは異なり、この誤差は系列長 \(T\) に依存しません。また、累積和・逆向き累積和・要素ごとの加算というテンソル演算だけで構成され、明示的な逐次ループもありません。qfeval-functions の msum・ma はこの方法で実装しています。

一方、各部分の和を単純に足すだけでは分散を求められません。分散は平均からの偏差の二乗和に基づく非線形な量であり、基準となる平均が窓ごとに異なるためです。次節では、各部分について保持する量を、和から分散の統計量へ拡張します。

2. 移動分散への拡張

分散には、2 本の移動和から求める方法がよく知られています。しかし、この方法は数値的に不安定です（§2.1）。一方、数値的に安定した既存の代替手法には、計算量が増える、または時点方向にベクトル化できないという問題があります（§2.2）。この節では、§1 のチャンク構造をそのまま利用し、各部分の和を分散統計量に置き換えます。これにより、移動和と同じく、計算量 \(O(T)\)・ベクトル化可能・誤差が局所的、という性質を持つ移動分散を構成します。

2.1. 二乗和の公式はデータの水準が大きいと桁落ちする

よく使われる \(O(T)\) の方法は、\(x\) と \(x^2\) の 2 本の移動和から分散を求める「二乗和の公式」です。

\[ \mathrm{Var} = \frac{1}{w - \mathrm{ddof}}\left(\sum_{i} x_i^2 - \frac{\bigl(\sum_{i} x_i\bigr)^2}{w}\right) \]

平均が \(\mu\)、分散が \(\sigma^2\) の窓を考えます。この式は、\(\sum x^2 \approx w(\mu^2 + \sigma^2)\) と \((\sum x)^2/w \approx w\mu^2\) という 2 つの大きな数の差から、小さな値 \(w\sigma^2\) を取り出します。両者の相対丸め誤差は \(\varepsilon\) のオーダーなので、結果の相対誤差はおよそ \(\varepsilon\,(1 + \mu^2/\sigma^2)\) に達します。\(\mu = 10^6\)、\(\sigma = 1\) のデータを float32 で計算すると、誤差は \(6 \times 10^4\) 倍に増幅され、有効桁が 1 桁も残りません。丸めによって差が負になることも多く、移動標準偏差を求めるために平方根を取ると NaN が生じます。実測でも、この見積もりどおりに精度が失われます（§3.3）。§1 の方法で移動和を安定に求めても、この減算による桁落ちは解消できません。原因は移動和の計算方法ではなく、二乗和の公式そのものにあります。

\(\sum x^2\) と \((\sum x)^2/w\) は、どちらも \(w(\mu^2 + \sigma^2)\) 程度の大きさです。浮動小数点数として表現するだけで、それぞれ \(\varepsilon\, w \mu^2\) 程度の絶対誤差を持ちます。一方、差の真値は \(w\sigma^2\) なので、相対誤差は \(\varepsilon\,\mu^2/\sigma^2\) のオーダーまで増大します。\(1 +\) の項は、\(\mu = 0\) の場合にも、丸めそのものに由来する相対誤差 \(\varepsilon\) が残ることを表します。

2.2. 安定な代替手法では速度かベクトル化が犠牲になる

各窓を独立に再計算すれば、数値的に安定した結果が得られます。例えば、平均を求めてから偏差二乗和を求める 2 パス法や、Welford 法が使えます。ただし、これらの計算量は \(O(Tw)\) です。PyTorch では、unfold で各窓のビューを作り、var を適用する形で実装できます。筆者の環境で測定したところ、\(O(T)\) の方法より 8〜300 倍遅くなりました（§4）。計算時間が窓幅に比例するため、\(w\) が数百になる用途には適しません。

もう 1 つの方法は、窓に入る要素を加え、窓から出る要素を取り除く逐次更新（Welford 法の add/remove 版）です。計算量は \(O(T)\) で、数値的にも比較的安定しています。pandas の rolling はこの方法を使っています。ただし、時点ごとに前の状態を引き継ぐため、時点方向には並列化できません。本記事では数千系列のバッチを float32 で処理することを想定しているため、時点方向にも並列化できるテンソル演算を採用します。

2.3. 部分ごとの統計量と局所中心化

ここで、チャンク分割（§1.2）に戻ります。分散を求めるため、各部分について要素数 \(n\)・平均 \(\mu\)・偏差二乗和 \(M = \sum (x - \mu)^2\) の 3 つ組を保持します。接頭辞の \(M\) を累積和から求めるには二乗和の公式が必要ですが、そのまま適用すると §2.1 の桁落ちが再発します。そこで、チャンク内の各値から基準値を引いた偏差に対して公式を適用します。基準値を \(r\)、偏差を \(y_i = x_i - r\) とします。分散と \(M\) は、値全体を平行移動しても変わりません。このため、任意の基準値 \(r\) について

\[ S = \sum_i y_i,\qquad Q = \sum_i y_i^2,\qquad \bar m = \frac{S}{n},\qquad M = Q - \frac{S^2}{n},\qquad \mu = r + \bar m \]

が成り立ちます。\(M = Q - S^2/n\) は §2.1 と同じ形の式ですが、\(y\) は基準値からの局所的な偏差です。そのため、水準 \(\mu\) によって桁落ちが際限なく大きくなることを防げます（§3.1）。接頭辞について \(y\) と \(y^2\) の累積和を、接尾辞について逆向きの累積和を計算します。この合計 4 本の累積和から、すべての部分の \((n, \bar m, M)\) を求められます。

基準値には、接頭辞ではそのチャンクの先頭の値を、接尾辞では末尾の値を使います。この選び方には、次の 2 つの利点があります。

誤差を窓内のスケールに抑えられる：どの接頭辞も先頭要素を含み、どの接尾辞も末尾要素を含むため、基準値は必ずその部分自身の要素です（図 2）。各部分は窓の部分集合なので、偏差 \(y\) の大きさは窓内の値域を超えません。系列の水準が高い場合やドリフトがある場合でも、計算に現れる値は窓の近傍における差分だけです（図 3・§3.1）。
NaN・±inf の影響を正確に限定できる：基準値が NaN や無限大の場合に影響を受けるのは、その要素自身を含む部分、つまり本来 NaN になるべき窓だけです（§3.2）。基準値にチャンクの平均のような集約値を使うと、チャンク内の 1 点にある NaN が、それを含まない窓にまで影響します。

基準値の選び方（w = 4）。隣接する 2 つのチャンクについて、長さ 1〜3 のすべての接尾辞と長さ 1〜4 のすべての接頭辞を並べた図。どの接尾辞もチャンク末尾の値 a を含み、どの接頭辞もチャンク先頭の値 b を含む。a と b は系列上で隣接している。 — 図 2：基準値の選び方（\(w = 4\)）どの接尾辞もチャンク末尾の値 \(a\) を含み、どの接頭辞もチャンク先頭の値 \(b\) を含むため、基準値は必ずその部分自身の要素になる。2 つの基準値 \(a\)・\(b\) は系列上で隣接し、その差 \(b - a\) は隣接要素の差になる（§2.4）。

局所中心化の前後でのスケールの比較。左は生の値のスケールで、値の列は 10 の 6 乗付近に平らに並び、先頭からそこまでの標準偏差は 0 に重なって区別できない。右は基準値 b を引いた後で、偏差は ±2 程度、そこまでの標準偏差は 1 程度と、どちらも窓内のスケールで見える。 — 図 3：局所中心化の前後でのスケール（水準 \(10^6\) の例）チャンク \(c\) の要素 \(x_5, \dots, x_8\)（オレンジ色）と、先頭からそこまでの標準偏差（灰色）を、生の値（左）と、基準値 \(b = x_5\) を引いた偏差（右）で示す。生の値のスケールでは、求めたい標準偏差（≈ 1）は 0 と区別できず、有効桁が約 7 桁の float32 では二乗和の差から取り出せない。基準値を引くと、偏差もそこまでの標準偏差も同じ窓内のスケールに現れる。標準偏差の値自体は、平行移動で変わらない。

2.4. Chan らのマージ式で 2 つの部分を結合する

2 つの部分 \(A\)（接尾辞）と \(B\)（接頭辞）の統計量を、冒頭に示した Chan らのマージ式[2]で結合します。

\[ M_{A\cup B} = M_A + M_B + \frac{n_A\,n_B}{n_A + n_B}\,\delta^2,\qquad \delta = \mu_B - \mu_A \]

この式は近似式ではなく恒等式です。\(n_B = 1\) とした場合は、Welford のオンライン更新式[1]と一致します。全体の平均を \(\mu = (n_A \mu_A + n_B \mu_B)/(n_A + n_B)\) とすると、各部分の寄与は「部分内のばらつき」と「部分間の平均の差」に分解できます。したがって、式の 3 項は順に \(A\) 内のばらつき、\(B\) 内のばらつき、\(A\) と \(B\) の間のばらつきを表します。

\(n = n_A + n_B\) とおき、\(x - \mu = (x - \mu_A) + (\mu_A - \mu)\) を二乗して \(A\) 内で足すと、\(\sum_{x \in A}(x - \mu_A) = 0\) によりクロス項が消え、\(\sum_{x \in A}(x - \mu)^2 = M_A + n_A(\mu_A - \mu)^2\) となります。\(B\) についても同様です。さらに \(\mu_A - \mu = -n_B\delta/n\)、\(\mu_B - \mu = n_A\delta/n\) を代入すると、追加項は \(\delta^2(n_A n_B^2 + n_B n_A^2)/n^2 = \delta^2 n_A n_B/n\) となり、マージ式が得られます。

\(\delta\) も局所的な値だけから計算できます。両方の部分の平均は「基準値＋相対平均」と表せるため、

\[ \delta = (b - a) + \bigl(\bar m_B - \bar m_A\bigr) \]

と書けます。第 1 項の \(b - a\) は 2 つの基準値の差で、\(a\) はチャンク \(c-1\) の末尾の値、\(b\) はチャンク \(c\) の先頭の値です。つまり、これらは系列上で隣接する 2 要素です。隣接要素の差を取ることでデータの水準は相殺され、その後の計算もすべて局所的な差分だけで行われます。接尾辞が空になる \(k = w\) の窓では、チャンク全体を接頭辞とする \(M_B\) をそのまま使います。最後に

\[ \mathrm{Var}_t = \frac{M_{A \cup B}}{w - \mathrm{ddof}} \]

として分散を求めます。標本分散なら \(\mathrm{ddof} = 1\)、母分散なら \(0\) です。移動標準偏差は、この分散の平方根です。

2.5. アルゴリズム全体

系列の先頭に NaN を追加して長さを \(w\) の倍数に揃え、長さ \(w\) のチャンクに区切ります（§1.2）。
各チャンクで、先頭の値を基準にした偏差 \(y\) について \(y\) と \(y^2\) の累積和を計算し、すべての接頭辞の \((n, \bar m, M)\) を求めます。
同様に、末尾の値を基準にした偏差について逆向きの累積和を計算し、すべての接尾辞の \((n, \bar m, M)\) を求めます。
各窓について、マージ式から \(M_{A \cup B}\) を求め、\(w - \mathrm{ddof}\) で割ります。窓と各部分の対応は添字をずらすだけで得られるため、この段階も要素ごとの演算で処理できます。
先頭に追加した NaN により、先頭の \(w-1\) 時点が自動的に NaN になります。\(T < w\) の場合は、すべての時点が NaN になります。

時点方向に依存関係がある演算は 4 本の累積和だけです。累積和は、GPU 上でも並列スキャンとして実装できる基本演算です。PyTorch では、累積和・逆向き累積和と要素ごとの演算だけを使い、30 行程度で実装できます（qfeval-functions の mvar）。

累積和は一見すると逐次計算に見えますが、部分和を木構造で組み合わせることにより、各チャンクの計算を深さ \(O(\log w)\) の並列計算に分解できます。この処理を prefix scan と呼び、並列ハードウェア向けの cumsum で広く使われています。

2.6. 小さな数値例

\(x = (50, 48, 49, 51, 52)\)、\(w = 3\) とします。先頭に NaN（以下では ∅ と表記します）を 4 個追加して長さを 9 にし、3 つのチャンクに区切ります。

\[ [\,∅,\ ∅,\ ∅\,]\quad[\,∅,\ 50,\ 48\,]\quad[\,49,\ 51,\ 52\,] \]

\(x_4 = 51\) で終わる窓 \((48, 49, 51)\) は、チャンク 2 の接尾辞 \(A = (48)\) とチャンク 3 の接頭辞 \(B = (49, 51)\) に分かれます。

部分	基準値	偏差 \(y\)	\(n\)	\(\bar m\)	\(M\)	\(\mu\)
\(A = (48)\)	末尾の値 48	\((0)\)	1	0	0	48
\(B = (49, 51)\)	先頭の値 49	\((0, 2)\)	2	1	\(4 - 2^2/2 = 2\)	50

基準値の差は隣接要素の差 \(49 - 48 = 1\) なので \(\delta = 1 + (1 - 0) = 2\) となり、マージ式から

\[ M_{A \cup B} = 0 + 2 + \frac{1 \cdot 2}{3} \cdot 2^2 = \frac{14}{3},\qquad \mathrm{Var} = \frac{14/3}{3 - 1} = \frac{7}{3} \]

を得ます。平均 \(148/3\) から直接計算した偏差二乗和も、\(16/9 + 1/9 + 25/9 = 14/3\) となり、一致します。系列全体に \(10^6\) を加えても、入力値を浮動小数点数へ丸める際の量子化を除けば、この計算に現れる値は基準値の差 \(1\) を含めて変わりません。これが局所中心化の効果です。

3. 数値安定性の性質

3.1. 誤差が窓の局所スケールに収まる理由

桁落ちが起こりうるのは、各部分で \(M = Q - S^2/n\) を計算する箇所です。ただし、基準値はその部分自身の要素なので、偏差 \(y\) の大きさは窓内の値域 \(R\) を超えません。値域が \(R\) である \(w\) 点の母分散には、\(\sigma^2 \ge R^2/(2w)\) という下界があります。これは、値域の両端にある 2 点だけでも偏差二乗和に \(R^2/2\) 以上寄与するためです。したがって、減算によって失われる桁数には上限があります。相対誤差は最悪でも \(O(\varepsilon w^2)\) であり、通常はそれより大幅に小さくなります。また、系列全体の水準 \(\mu\)・ドリフト・系列長 \(T\) には依存しません。これは、二乗和の公式の誤差因子 \(\mu^2/\sigma^2\)（§2.1）が水準とともに際限なく増大することと対照的です。

\(|y| \le R\) なので \(Q = \sum y^2 \le wR^2\) です。\(S\)・\(Q\) の累積和と後続演算の絶対誤差は最悪でも \(O(\varepsilon w^2 R^2)\) ですが、真値は \(M = w\sigma^2 \ge R^2/2\) なので、相対誤差は \(O(\varepsilon w^2)\) で抑えられます。実測（§3.3）では float32・\(w = 20\)〜\(50\) の条件でアルゴリズム由来の誤差は最大 \(10^{-6}\) 未満でした。窓ごとの 2 パス法や Welford 法の誤差は \(O(\varepsilon w)\) であり、極端に大きな窓幅で最後の 1 桁まで詰めたい場合は 2 パス法が優位です（§7）。

また、窓内の値がすべて等しい場合、偏差 \(y\) は厳密に 0 となるため、分散も丸め誤差なしで厳密に 0 となります。二乗和の公式では、値が一定の窓でも丸め誤差が生じ、負の値になる場合があります。

3.2. NaN・±inf の影響を窓内に限定できる理由

「窓に NaN や無限大が含まれていれば結果は NaN になり、含まれていなければ影響を受けない」という性質が、厳密に成り立ちます。NaN や無限大を含まない窓の結果は、それらを別の有限値に置き換えた場合とビット単位で一致します。この性質は、基準値として各部分自身の要素を選ぶことで得られます。実装のテストでは、NaN と ±inf を各時点に順に置き、すべての位置でこの性質が成り立つことを検証しています（§6）。

NaN や無限大が影響する経路は 3 つあります。（1）各部分の累積和：ある位置までの累積和は、それより後の要素に依存しません。そのため、該当する値より手前で終わる接頭辞と、該当する値より後ろから始まる接尾辞の統計量は影響を受けません。（2）基準値：接頭辞の基準値であるチャンク先頭の値は、すべての接頭辞に含まれます。同様に、接尾辞の基準値であるチャンク末尾の値は、すべての接尾辞に含まれます。したがって、基準値が NaN や無限大であれば、その部分は必ず該当する値自体を含んでいます。（3）マージ式の \(\delta\) が参照する 2 つの基準値も、それぞれ接尾辞と接頭辞に含まれる要素です。以上から、影響が及ぶ範囲は「NaN や無限大を要素として含む窓」と正確に一致します。

3.3. 実測

float32 を使い、本手法と二乗和の公式（§2.1）の精度を比較しました。参照値には、同じ系列を float64 で計算した結果を使います。「オフセット」は平均 \(10^6\)、分散 1 の乱数系列（\(w = 20\)）です。「ドリフト」は、水準 1000 の周辺を標準偏差 0.01 の小さなステップで変動するランダムウォーク（\(w = 50\)）です。後者は、大域的な定数を 1 つ引くだけでは精度を改善できない例です。

ケース（float32）	二乗和の公式	本手法
オフセット：相対誤差（中央値）	約 \(9 \times 10^{4}\)（有効桁が残らない）	\(6 \times 10^{-3}\)（入力の量子化誤差が下限）
オフセット：負の分散の割合	31%	0
ドリフト：相対誤差（中央値）	約 \(72\)	\(1.3 \times 10^{-4}\)

本手法で残る誤差は、主に float32 で入力を表現するときの量子化によるものです。\(10^6\) 付近における float32 の刻み幅は約 0.06 です。同一の float32 入力を float64 で再計算した参照値と比べると、アルゴリズム自体に由来する誤差は、どちらの場合も最大 \(10^{-6}\) 未満でした。

4. 速度

二乗和の公式が 2 本の累積和と要素ごとの演算で構成されるのに対し、本手法では累積和が 4 本に増え、要素ごとの演算も多くなります。一方、窓ごとの再計算（unfold + var）は計算量が \(O(Tw)\) であり、窓幅に比例して実行時間が増えます。筆者の環境（Apple シリコンの CPU、PyTorch 2.7、float32）で測定した結果は次の通りです。

形状	窓幅	二乗和の公式	窓ごと再計算	本手法
(10000,)	20	0.21 ms	1.5 ms	0.18 ms
(100, 10000)	20	2.0 ms	145 ms	2.7 ms
(100, 10000)	500	2.1 ms	763 ms	3.1 ms
(1000, 10000)	500	21 ms	7237 ms	23 ms

本手法の実行時間は二乗和の公式の 0.9〜1.5 倍に収まり、窓幅を 20 から 500 に広げてもほとんど変わりません。窓ごとの再計算と比べると、上の表の範囲では、1 系列の場合に 8 倍、バッチ処理の場合に 50〜300 倍高速です。GPU でも効率よく実行できる構成ですが、表に示したのは CPU での結果であり、GPU 上の速度は検証していません。

5. 関連手法との比較

Welford のオンラインアルゴリズム[1]は、要素を 1 つずつ追加するときに平均と分散統計量を安定して更新する方法です。Chan らの式[2]は、これを任意の 2 つの部分集合を結合できるように一般化したものです。この式は、もともと並列計算やペアワイズ集計のために使われてきました。本手法では、「どの移動窓も常に 2 つ以下の部分に分けられるチャンク分割」とこの式を組み合わせ、すべての時点を一括のテンソル演算で処理します。また、総和の丸め誤差をブロック分割によって抑えるペアワイズ加算の考え方を、移動窓に合わせて 1 段だけ適用したものとも解釈できます。各方式の違いは次の通りです。

ペアワイズ加算は、要素列を前半と後半に分け、それぞれの和を再帰的に計算する方法です。逐次加算では丸め誤差が項数 \(T\) に比例して蓄積する可能性がありますが、ペアワイズ加算では、加算を表す木の深さ \(\log T\) に比例する程度に抑えられます。numpy の sum などで採用されている標準的な手法です。

方式	計算量	ベクトル化	数値安定性
二乗和の公式（移動和 2 本）	\(O(T)\)	可	不安定（誤差因子 \(\mu^2/\sigma^2\)）
大域累積和の差分	\(O(T)\)	可	系列長に比例して劣化（§1.1）
窓ごとの再計算	\(O(Tw)\)	可	安定（\(O(\varepsilon w)\)）
逐次 add/remove（pandas など）	\(O(T)\)	時点方向は不可	概ね安定（除去の実装に依存）
本手法（チャンク＋局所中心化＋マージ）	\(O(T)\)	可	安定（最悪 \(O(\varepsilon w^2)\)、水準・系列長に非依存）

分散を数値的に安定して計算する方法については、各公式で精度が低下する条件や、データに近い値によるシフトを含め、Chan らの分析[3]と Higham の教科書[4]に整理されています。本手法の局所中心化は、そこで推奨されているシフトを、各部分について自動的に選ぶ方法です。この選び方では、NaN の影響範囲も正確に保たれます（§2.3）。

6. 実装メモ

先頭に追加する NaN の役割：チャンク境界を揃えるために系列の先頭へ追加する値を NaN にすると、「先頭の \(w-1\) 時点は NaN」「\(T < w\) なら全時点が NaN」という出力規約を、専用の分岐なしで満たせます。NaN を \(w\) 個以上追加することで、実データの先頭に対しても直前のチャンクを用意できます。
ddof の扱い：\(0 \le \mathrm{ddof} < w\) を前提とします。通常は標本分散の \(1\) か母分散の \(0\) です。
移動平均・移動標準偏差：移動和と移動平均は §1 のチャンク和だけで求められます。移動標準偏差は、移動分散の平方根として求めます。
負値のクランプ：本手法の \(M\) は、丸め誤差によってごくわずかに負になる可能性があります。一方、通常の 2 パス法では非負であることが保証されます。平方根を取る場合など、厳密な非負性が必要であれば、最後に 0 を下限としてクランプします。実測の条件（§3.3）では、負の値は観測されませんでした。

実装の正しさを継続的に確認するため、次の 5 項目を回帰テストに含めました。（1）系列長と窓幅のすべての組み合わせについて、チャンク境界の位置にかかわらず、窓ごとの再計算の結果と一致すること。（2）大きなオフセットを持つ float32 の系列で、同一入力を float64 で計算した結果との相対誤差がしきい値以下であること（§3.3）。（3）ランダムウォーク系列でも同じ条件を満たすこと。（4）NaN と ±inf を各時点に順に置き、その影響が該当する値を含む窓だけに限定されること（§3.2）。（5）値が一定の窓では厳密に 0 を返し、dtype が保たれること。

7. 設計上の限界

誤差の上界が 2 パス法より大きい：最悪の場合の誤差は \(O(\varepsilon w^2)\) であり、窓ごとの 2 パス法における \(O(\varepsilon w)\) より \(w\) 倍大きい上界です。速度よりも末尾の桁までの精度を優先する場合や、\(M\) が必ず非負である必要がある場合は、2 パス法が適しています（§3.1・§6）。
ストリーミング処理には適さない：本手法はチャンク単位の一括計算を前提としています。1 点ずつ到着するデータを低遅延で逐次処理する場合は、Welford 法や add/remove 方式（§2.2）の方が適しています。本手法はバッチ計算向けです。
NaN を無視して計算できない：窓に NaN が含まれていれば、その窓の結果も NaN になります。pandas の min_periods のように、NaN を除いた要素だけで計算する必要がある場合は、有効な要素数のカウントを含む別の構成が必要です。
入力は浮動小数点数に限られる：先頭への NaN の追加と偏差の計算を行うため、整数系列はあらかじめ浮動小数点数へ変換する必要があります。

8. まとめ

大域累積和の差分による移動和は、系列が長くなるほど精度が低下します。系列を窓幅と同じ長さのチャンクに区切ると、どの窓も「接尾辞＋接頭辞」という最大 2 つの部分に分けられます。各部分の和は最大 \(w\) 項なので、移動和と移動平均を安定して求められます。
移動分散を二乗和の公式で求めると、誤差が \(\mu^2/\sigma^2\) に比例して増幅されます。一方、安定した窓ごとの再計算には \(O(Tw)\) の計算量が必要であり、逐次更新は時点方向に並列化できません。本手法では、各部分について保持する量を、和から「要素数・平均・偏差二乗和」に置き換えます。統計量をその部分自身の要素を基準にした偏差で表すことで、計算全体を局所的な差分だけで構成し、Chan らのマージ式で 2 つの部分を結合します。
これにより、窓幅に依存しない計算量 \(O(T)\)、テンソル演算だけによる実装、水準・ドリフト・系列長に依存しない局所的な誤差、NaN と無限大の正確な伝播を同時に実現できます。実測した処理時間は、二乗和の公式の 0.9〜1.5 倍でした。
用いる要素技術は、1960〜80 年代に発表された Welford と Chan–Golub–LeVeque の古典的な手法です。本記事では、これらを移動窓のテンソル計算に組み合わせる構成を整理しました。実装は qfeval-functions の msum・mvar で公開しています。

9. 参考文献

オンライン更新式・マージ式の原典と、分散計算の数値誤差解析の標準的な文献は以下の通りです。

[1] B. P. Welford, "Note on a Method for Calculating Corrected Sums of Squares and Products," Technometrics 4(3), 419―420, 1962. （Welford のオンライン更新式。逐次追加に対する平均と偏差二乗和の安定な更新）
[2] Tony F. Chan, Gene H. Golub, Randall J. LeVeque, "Updating Formulae and a Pairwise Algorithm for Computing Sample Variances," COMPSTAT 1982, 30―41, 1982. （Chan–Golub–LeVeque による更新式とペアワイズアルゴリズム。本記事のマージ式の原典）
[3] Tony F. Chan, Gene H. Golub, Randall J. LeVeque, "Algorithms for Computing the Sample Variance: Analysis and Recommendations," The American Statistician 37(3), 242―247, 1983. （同著者らによる分散計算アルゴリズムの誤差分析と推奨。シフト（中心化）の効果を含む）
[4] Nicholas J. Higham, "Accuracy and Stability of Numerical Algorithms, Second Edition," Society for Industrial and Applied Mathematics, 2002. （丸め誤差解析の標準的教科書。総和・分散計算の誤差上界）

最終更新日: 2026年7月16日