ストリームアルゴリズム
この記事はレイトレ合宿5!?のアドベントカレンダーの7週目の記事です。
無限の大きさがあるデータの集合から、小さなフットプリントで要約統計量を得たい。
ストリームアルゴリズム
何らかの集まりに対して、そこからの何らかの統計量を得たいと思うことがよくあります。例えば数値の列から平均を求めるには？と言われたら次のようなコードを書くかもしれません
では今度は数が決まっていないデータが来た場合はどうでしょうか。この例で出てきた配列は用意できません。しかし少し考えると次のようにすれば算出できることに気が付きます。
これがストリームアルゴリズムです。全てをメモリ上に収めるのが難しい集合でも、このようにある種の要約統計量を求めることができる場合があります。ストリームとしてデータがやってくるようなソフトウェアにはどのようなものがあるでしょうか？よく考えるとたくさんあります。プレイヤーからコントローラー入力を受け続けるゲーム、HTTPレスポンスを処理するブラウザ、センサー値を発し続けるIoT的なもの、そしてレイトレです！
レイトレを実装するうえで便利なストリームアルゴリズムをいくつか紹介します。
最大値最小値
最初はとても簡単な例、最大値、最小値です。最大、最小の操作は結合法則が成り立ちます。
$$ max(a,b,c) = max(max(a,b), c) $$
つまり、「それまでの最大値と次の値との最大値」は「集合全体の最大値」になります。コードにすると次のようになります。
当たり前ですね。
使用例
あまりにも無意識に使っていることが多いと思われますが、いろんなところに使っているかと思います。Fireflyが出ていないかをチェックするために最大輝度を出したり、ロシアンルーレットがバグってやたら長いパスが出来ていないかをチェックする…などなど。
相加平均
次は平均(相加平均)を出してみます。相加平均$\mu$の定義は次の通りです。
$$ \mu = \frac{1}{n} \sum_{i=1}^{n} {x_i} $$
ここで$i$番目のサンプルが$x_i$で、全てのサンプルの数は$n$です。これは最初の例のように、それまでの総和を蓄えておいて、平均が必要になったらそれまでの要素数で割れば良さそうですし、少し工夫するとsumを消すこともできます。
使用例
「平均」は最も日常で触れることが多い統計量で、よく慣れ親しんでいるので、大量の数値が流れていくところのバグを追う最初の手がかりに使ったり、またある時はサブピクセルの輝度の平均(box filter)のように直接使う場合もあると思います。
分散、標準偏差
分散もストリームアルゴリズムで求めることができます。確率変数を$X$とするとその分散$Var(X)$の定義は次の通りです。
$$ Var(X) = E[(X-E[X])^2] $$
ここで$E[]$はその中の確率変数の期待値を表します。$X$の期待値(つまり平均値)からの偏差の二乗の期待値です。具体的なサンプルからなる標本分散$s^2$は、$i$番目のサンプルを$x_i$とすると次のようになります。
$$ s^2 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \mu)^2 $$
この標本分散は母分散の一致推定量ではあるのですが、不偏推定量ではありません。標本分散に$\cfrac{n}{n-1}$を掛けると母分散の不偏推定量の不偏分散になり、こちらもよく使われます。
$$ \hat{\sigma}^2 = \frac{1}{n-1}\sum_{i=1}^{n} (x_i - \mu)^2 $$
さてこの標本分散と不偏分散ですが、$\mu$はその要素数までの相加平均なので、先に$\mu$が分かっていないと、つまり最初に全てのサンプルに対する相加平均が分かっていないと出しようにないように一見、見えます。疑似コードは次のようになります(これは動きません)。
しかしいくつかの式変形を経ると驚くことに逐次的に求めることができることがわかります。詳しくは"分散の逐次更新“を参照してください。次のコードは不偏分散を逐次的に出しています。
使用例
不偏分散はかなり強力な値でここからまず母平均の区間推定が出来ます。十分な数のサンプルがあるとき、母平均が次の区間の間に収まる分布である確率が95%になります。(95%の信頼区間)
$$ \mu0 - 1.96 \sqrt{\cfrac{\hat{\sigma}^2}{n}} \leq \mu \leq \mu0 + 1.96 \sqrt{\cfrac{\hat{\sigma}^2}{n}} $$ ※ 1.96は自由度無限のt分布の両側5パーセント点より。
CLTが成立することからもわかるように、サンプル数が十分に多ければ、任意の母分布でこれは成立します。母平均の信頼区間が出ると、レイトレでのピクセルの値がどれほど収束しているかがわかります。もし母平均の信頼区間が出力先の精度のオーダーを下回る、例えば8bit深度であれば信頼区間が1/255を下回り始めたらそれ以上サンプルしても結果は大きく変わらないかもしれません。逆に信頼区間が広く全く収束していないと予想されるピクセルがあればそこに計算資源を多めに割くことでイメージ全体のクオリティを上げることができるかもしれません。
また別の使い道としては、外れ値の判定があります。母分散を$\sigma^2$とするとチェビシェフの不等式より分布に従うあるサンプルが存在する確率の上限が出ます。
$$ Prob(|X-\mu| \geq k \sigma) \leq \cfrac{1}{k^2} $$
「偏差値70以上と30以下の人数の和は全体の1/4を超える人数には絶対にならない」と読み替えるとわかりやすいかもしれません。かなり緩い条件であることからもわかるように、これも任意の母分布で成立します。これでとても低い確率でしか出えないサンプル、レイトレで言えば超高輝度なサンプルですが、これが判定ができるようになり、それに対して何らかの特殊な操作をすることで輝点を防ぐこともできます。このようなことをするとレンダリング結果にバイアスが乗ってしまいますが、実際に分散をクランプしてしまうようなレンダラーもあります。
共分散、ピアソンの積率相関係数
分散と同様に共分散も計算できます。共分散の定義は次の通りです。
$$ Cov(X,Y) = E[(X-E[X])(Y-E[Y])] $$
具体的なサンプルからなる標本共分散は次のようになります。
$$ s = \cfrac{1}{n} \sum_{i=0}^{n} ({x_i - E(X)}) ({y_i - E(Y)}) $$
これもまた分散の時と同じように逐次的に求めることができます。
さて共分散と分散が分かれば、ピアソンの積率相関係数がわかります。よく「相関係数」と呼ばれる値です。共分散についていた次元がとれて無次元量になっており-1から1の間を取り、相関している度合いをある程度は直観的に説明、理解できます。ピアソンの積率相関係数の定義は次の通りです。
$$ \rho = \cfrac{Cov(X,Y)}{\sqrt{Var(X)Var(Y)}} $$
使用例
相関係数が分かるので、二つの確率変数が正規分布に従えばt検定が、そうかわからなければPearsonの無相関検定を使って、二つの確率変数が相関するとは言えないかがわかります。サンプラーが壊れていないか(異なる次元が相関してしまっていないか)をチェックするのに使えるのではないでしょうか。
最頻値
最大、最小、平均、分散、共分散が求まったので、次は最頻値です。しかし最頻値を固定のメモリフットプリントで確実に得る方法は存在しません。例えば100億個のサンプルがあり、一つの重複のみを除き全て種類が違うサンプルの場合、99億9999万9999個の頻度を保存する領域が必要になることからも理解できるかと思います。
ただし事前に決めたフットプリントで決定する誤差を許容することで、その誤差以内の頻度を保証するLossy Countingという方法があります。
- まずストリームデータを固定サイズのバケットに分割して考えます。
- 頻度を保存するテーブルを用意します。そこにIDと度数と頻度を最後に更新したバケットの番号を記録します。
- 一つのバケットを処理するたびに十分小さな頻度の要素をテーブルから削除します。「小さな頻度」の基準は、度数とその要素を最後に更新したバケット番号の和が、現在のバケット番号+1よりも小さいかです。
- テーブル内にある要素は常に、頻度-バケット番号 以上の頻度が出てきたことが保証されます。
頻度を保存するテーブルも常にバケットのサイズを超えることがありません。コードにすると次の通りです。
使用例
レイに当たったオブジェクトのIDなどの頻度を出せば、意図通りのオブジェクトに当たっているか(当たっていないか)をチェックできます。複数の戦略をとるようなアルゴリズムの場合にどれがどれほどの割合で採用されたかなどもわかるかもしれません。
Reservoir sampling
最後に少し毛色が違う手法です。特定のパラメーターの統計量を得るのではなく、無限の大きさの集合から、その集合を代表するような小さな集合を作り出すことが出来るReservoir samplingの紹介です。
- まずサンプルを入れる配列を用意します。
- 配列の要素数までのデータは前から順番に格納していきます。
- 配列の要素数を超えるデータが来たら、[0,データ数)の範囲の乱数を作成し、それをインデックスとして配列に格納します。インデックスが配列の範囲外になったら棄却します。
使用例
目grepをしたいが、膨大過ぎて全てのデータに目を通せず、またデータが場所によって偏りもあるような場合でも、全体像をつかみたい場合など。
最後に
どんなに大きな集合でもその母集団から情報を取り出すいくつかの方法の紹介をしました。ハマるとつらいレイトレのデバッグもこれらの道具を使って楽になればいいですね。