この記事はレイトレアドベントカレンダー(2018)の記事です。

はじめに

コンピュータープログラムには、「めったにないけれど、万が一が起きた時にする処理」(フォールバック)が用意されていることがあります。例えば、レイトレーシングでは、「あまりにもグレージングアングルすぎるからなかったことにする」のような処理です。

float pdf(Vec3 wi)
{
	// wiは入射方向でzは[0,1)であることを期待する
	// 小さすぎるzは巨大な値を返してしまうので握りつぶす
	if(wi.z < 1e-8f)
	{
		return 0.0f;
	}
	else
	{
		return 1.0f / wi.z;
	}
}

このような処理はめったに起こらないがために最終結果に影響を与えないという場合が多く、「めったに起こらない」という仮定が崩れると違った結論が導かれてしまうことがあります。例えばこの例のwiがZが負の領域も含めて返していたらこの関数の意図通りの処理となるでしょうか? このように稀であってほしいイベントが本当に稀なのかを確認する手段が必要です。

ここからは記号がいくつか出てくるので先に整理をします。

お役立ち情報コーナー サンプル数 vs サンプルサイズ

CGの文脈では「観測値が3個手に入った」ことを「サンプル数3」と呼び、レンダラーの設定値にもそう書かれていることが多いですが、統計の文脈では「サンプルサイズ3」と呼び、「サンプル数」には「観測値が集まった集合の数」という別の意味があります。このことについてVeachも言及しています。

Note that there are some differences in the standard terminology for computer graphics, as compared to statistics. In statistics, the value of each $X_i$ is called an observation, the vector ($X_1$…$X_N$ ) is called the sample, and $N$ is called the sample size. In computer graphics, on the other hand, typically each of the individual $X_i$ is referred to as a sample, and N is the number of samples. We will normally use the graphics conventions.

つまりサンプル(標本)に「観測値がいくつか入ったもの」の意味がすでにあるので、レンダラーの"SPP(Sample Per Pixel)“という設定値名は統計学的には誤用です。しかしCGではCGの文脈に従い、観測値=サンプルとします。

ベルヌーイ試行→ベルヌーイ分布→二項分布

何らかのイベントが起こった場合を**“head”**、起こらなかった場合を**“tail”**と呼んでみます(コインの表裏をこう呼称します)。このように結果がhead/tailの二通りしかない試行のことを**ベルヌーイ試行(Bernoulli trial)**と呼び、この試行の結果の確率分布を**ベルヌーイ分布(Bernoulli distribution)**といいます。headになる確率がpだとすると次のように定義されます。

$$\begin{eqnarray} P[X]=\left \{ \begin{array} \\
p & if \ head \\
0 & if \ tail \end{array} \right. \end{eqnarray}$$

headを1、tailを0とするとこの確率分布の期待値は次のようになります。

$$\begin{eqnarray} E[X] &=& 1 \times P(head) + 0 \times P(tail) \\
&=& 1 \times p \\
&=& p \tag{1}\label{1} \end{eqnarray}$$

分散は次のようになります。

$$\begin{eqnarray} Var[X] &=& E[X^2] - E[X]^2 \\
&=& (1^2 \times P(head) + 0^2 \times P(tail)) - (p)^2 \\
&=& p(1-p) \tag{2}\label{2} \end{eqnarray}$$

ベルヌーイ分布はベルヌーイ試行を一回行った結果ですが、このベルヌーイ試行を何回も行った結果の**「headの数の総数の分布」**にも名前がついており、**二項分布(Binomial distribution)**と呼ばれています。二項分布は次のような確率度数関数になります。

$$\begin{eqnarray} P[X=k] &=& \binom{n}{k} p^k (1-p)^{n-k} \tag{3}\label{3} \end{eqnarray}$$

$n$が試行回数、$k$がheadであった回数です。少しイカツイ形をしていますが、これを直接扱わなくても、ベルヌーイ試行はお互いに独立した試行なので、期待値と分散は単にベルヌーイ分布の期待値$(\ref{1})$と分散$(\ref{2})$に試行回数nをかけたものになります。楽ですね。

$$\begin{eqnarray} E[X] &=& np \tag{4}\label{4} \\
Var[X] &=& np(1-p) \tag{5}\label{5} \\
\end{eqnarray}$$

二項分布→正規分布

二項分布の確率度数関数を使えば、各度数の確率が直接出せるので、仮定した頻度(「このif文は多くても0.01%未満の確率で起こる」など)よりも多く起こっているのかを計算すればよいのです。しかし一つ問題があります。$(\ref{3})$式中の二項係数の部分を効率的に実装できないのです。巨大な2次元のテーブルを作れば回避できるかもしれませんが、もう少しスマートに解決したいところです。

ところで二項分布は独立した試行を繰り返したものの和から得られたものですが、何か聞き覚えがあります。中心極限定理です。中心極限定理(Central limit theorem,CLT)は、(ほぼ)任意の母分布において、独立したn個をサンプルした場合に、nが大きければ、サンプルの和が正規分布$N(n\mu,n\sigma^2)$に近づいていくというものです。ベルヌーイ分布の期待値は$p$で、分散は$p(1-p)$だったので、二項分布は$N(np,np(1-p))$に近づいていくということになります。

正規分布は巨大なテーブルを作らなくても計算ができるので、便利になりました。

二項分布→ポワソン分布→正規分布

さてベルヌーイ試行の繰り返しを正規分布として仮定する話を書きましたが、今回の場合はまた別の分布があります。二項分布はhead/tailの試行結果が、大きな粒度で感じられる世界観で語られていましたが、nが非常に多い、つまり試行回数が非常に多い状態を考えてみます。このように非常にサンプル数が多い世界では各々の試行は細かい粒のようになり、簡略化した分布関数に置き換えることができます。$k$は有限の数のままで、$n$を無限大だと考えた場合の二項分布を考えてみます。

$$$$\begin{eqnarray} \lim_{n \to \infty} P[X=k] &=& \lim_{n \to \infty} \binom{N}{k} p^k (1-p)^{n-k} \\
&=& \lim_{n \to \infty} \cfrac{n!}{k!(n-k)!} \biggl( \cfrac{\lambda}{n} \biggr) ^k \biggl(1- \cfrac{\lambda}{n} \biggr)^{n-k} (変数 \lambda = npを導入) \\
&=& \cfrac{\lambda^k}{k!} \lim_{n \to \infty} \cfrac{n}{n} \cfrac{n-1}{n}\cfrac{n-2}{n}\cdots \cfrac{n-k+1}{n} \biggl(1- \cfrac{\lambda}{n} \biggr)^{n-k} \\
&=& \cfrac{\lambda^k}{k!} \lim_{n \to \infty} \biggl(1- \cfrac{\lambda}{n} \biggr)^{n-k} \\
&=& \cfrac{\lambda^k}{k!} \lim_{n \to \infty} \biggl(1- \cfrac{\lambda}{n} \biggr)^{n} \\
&=& \cfrac{\lambda^k e^{-\lambda} }{k!} \tag{6}\label{6} \\
\end{eqnarray}$$$$

この分布はポワソン分布(Poisson distribution)と呼ばれる分布で、期待値と分散が共に$\lambda$(導出割愛)となる面白い分布です。二項分布の連続関数版と考えることができます。

※ よく「ポワソン分布は稀な事象に対して適用される」と説明されますが、「各々の事象が細かい粒のように見えるほどnが巨大な時の分布」と説明した方が個人的にはかなりしっくりします。「ポワソンの少数の法則」という名前から誤解が広まった？ 今回のようにpが一定の時は「稀な事象を扱う関数」として、npが一定の時は「連続時間的な関数」として、接するのがいいのではないかと思っています。

統計量z

正規分布で近似できたものは、母平均の検定ができます。母平均の検定をするために次のような統計量$z$を導入します。

$$ z = \cfrac{m - \mu}{\sqrt{n \sigma^2}} \tag{7}\label{7} \\
$$

この統計量$z$の値の絶対値が大きいほど、より想定よりも外れた値が多いということになります。分子の平均の差が大きいほど、分母の分散とサンプル数の積が小さいほど、外れ値の可能性が上がるということからも直感的にわかりやすいかと思います。今回だと、「想定した稀さよりもずっと少ない」か「想定した稀さよりもずっと多い」です。前者は今回は問題ならないので後者の「想定した稀さよりもずっと多い」にのみ注目します(片側検定)。

手順としては次のようになります。

1.帰無仮説と対立仮説を立てる。今回は帰無仮説は「想定した頻度未満でしかイベントは起こらない」、対立仮説は「想定した頻度以上にこのイベントは起こる」です。

2.有意水準を設定します。「有意水準5%」としたら、「5%の確率で誤った結論が導かれる」ということになります。適当に有意水準1%として、正規分布表から引っ張ってくると2.33となるので、Z値が-2.33を下回る数字が出れば、帰無仮説が棄却されます。ただ厳密に棄却する/しないを知りたいというよりも確実に想定よりも稀なイベントだということを確認するという目的があるのでZ値は適当に-3などにしてもいいかもしれません。

統計量$z$の式$(\ref{7})$には、期待値と分散が必要なので、ここまでできた分布のまとめを載せます。

二項分布の統計量$z$は次のようになります。

$$\begin{eqnarray} z &=& \cfrac{m - np}{\sqrt{n^2 p(1-p)}} \\
\end{eqnarray}$$

同様にポワソン分布の統計量$z$は次のようになります。

$$\begin{eqnarray} z &=& \cfrac{m - np}{\sqrt{n^2 p}} \\
&=& \cfrac{m - \lambda}{\sqrt{n \lambda}} \\
\end{eqnarray}$$

これらから次のことがわかります。

• pが小さければ、二項分布もポワソン分布も統計量$z$は変わらない。
• ポワソン分布はpがわからなくてもnp(全体で発生した件数)があれば使える。

まとめ

稀なイベントが本当に稀であるのかの確認には次のように場合分けをするといいかもしれません。

• サンプル数: 少
  　→ 二項分布を直接使う
• サンプル数: 多 かつ 頻度の上限が見当が付いている
  　→ 二項分布の正規分布近似を使って、Z検定(or t検定)を行う
• サンプル数: 多 かつ 全体で起こる回数の上限が見当が付いている
  　→ ポワソン分布の正規分布近似を使って、Z検定(or t検定)を行う

サンプルコード

稀なイベントが発生しすぎていることがちゃんと補足されるかを次のそれぞれの場合で比較します。

1. 二項分布
2. 二項分布を正規分布で近似
3. ポワソン分布の正規分布近似

TODO: 後でちゃんと書きます。

参考

[1] Matt Pharr. CHECK_RARE and making sense of unusual occurrences. https://pharr.org/matt/blog/2018/05/31/check-rare.html
[2] Eric Veach. “Robust Monte Carlo Methods for Light Transport Simulation”.