このブログポストはレイトレアドベントカレンダー2019の23日目の投稿です。
はじめに¶
MCレイトレーシングでは、pdf関数とsample関数がBSDFをはじめ様々なクラスで出てきます。
class BSDF
{
public:
struct Sample { vec3 wi; float pdf; };
float pdf(vec3 wo, vec3 wi) const;
Sample sample(vec3 wo, vec2 uv) const;
};
pdf関数とsample関数は別々の関数として実装されるので、計算が間違っていたり、バグらせてしまったりすると、一貫性がなくなってしまいます。例えば、pdf関数で高い値が出ている領域にもかかわらず、sample関数では全くサンプルされないといった場合です。pdf関数とsample関数が一貫性のある妥当な実装になっているのかを検証する手段についてこのブログポストでは述べていきます。
Sin関数¶
まずは簡単な例として、sin関数のpdf関数とsample関数を実装してみます。
import math
# sinのpdf関数とsample関数
s = math.pi/2
def sin_pdf(x): return np.sin(x*s) * s
def sin_sample(u): return np.arccos(1.0-u)/s
このSin関数のpdf関数とsample関数をプロットしてみます。pdf関数は折れ線グラフで、sample関数はヒストグラムでプロットします。
import matplotlib.pyplot as plt
import numpy as np
import random
np.random.seed(1234)
random.seed(1234)
# pdf関数とsample関数からプロットを作る
def plotSamples(pdf,sample,n):
u = np.random.random(n)
samples = sample(u)
plt.xlim(0.0,1.0)
plt.ylim(0.0,1.7)
plt.hist(samples, bins=50, density=True, color='cornflowerblue', edgecolor='black')
x2 = np.linspace(0.0,1.0,100)
y2 = pdf(x2)
plt.plot(x2,y2, c='red')
plt.show()
plotSamples(sin_pdf,sin_sample, 10000)
目視でみてみると、一致しているようにみえます。しかしその程度の確認で本当に妥当さを確認したと言えるでしょうか？ 次のプロットは、pdf関数をsin関数から二次関数にしたものです。sin関数は二次関数ではないので、正確に一致はしないはずですが、目視ではかなり一致しているように見えるのではないでしょうか。
# Sin関数っぽく見えるが、ただの二次関数なので、sin_sample関数と一致した結果にはならない
def sin_like_pdf(x): return -3/2*(x-2)*x
plotSamples(sin_like_pdf, sin_sample, 10000)
nを増やせばズレが目立つようになりますが、それでも目で判別できないズレは見逃してしまいます。もう少し機械的、論理的、統計的にpdf関数とsample関数が一貫性がある実装になっているかを判定する方法が欲しくなってきます。
適合度の$χ^2$検定¶
何かの観測値のサンプルが、その理論的な発生する確率と一致しているかの検定が適合度の$χ^2$検定です。$χ^2$検定には母分散の検定もありますが、それとは別物であることに注意してください。
$k$種類のカテゴリー$A_1,A_2,...A_n$の観測度数$f_i$が理論確率$p_i$と一致、つまり$f_i=n_i p_i$の関係がちゃんと成り立っているかを検定します。
次のように定義した理論値と観測値とのズレをピアソンの適合度基準と呼びます。
$$ \chi^2 = \sum_{i=1}^{k}{\frac{(f_i-np_i)^2}{np_i}} $$式を眺めると、分子に理論値と観測値の差を表し、分母でそれを標準化のようなことをしていることがわかります。この確率変数$\chi^2$は自由度が$k-1$の$\chi^2$分布に従います。言葉で説明するよりも、具体的なコードをみたほうが伝わるので計算してみます。scipyに便利な関数があるのでそれをそのまま使います。
from scipy import stats
# 観測される各カテゴリのサンプル数
observed = [468, 494, 532, 494, 517, 495]
# 理論上の各カテゴリのサンプル数。どのカテゴリも等しい数サンプルされることを期待する。
expected = [500, 500, 500, 500, 500, 500]
# scipy.statsを使って$\chi^2$検定をする
result = stats.chisquare(observed, expected)
print(f'χ2={result.statistic} P-Value ={result.pvalue}')
P値(P-Value)とは、ある確率変数が従うはずの分布の中でどのあたりにいるかを示す値です。これが極端な値を取ると「従うはずの分布」という前提が怪しいということになります。今回のピアソンの適合度基準は、$\chi^2$分布に従うので、$\chi^2$分布の中でどのあたりに確率変数$\chi^2$があるかを示します。とても低い値がでてくれば、奇跡的な偶然(なので恐らく理論上のモデルからは大きく乖離している)となります。
このサンプルコードではP値は0.43と高い値が出ているので、期待した分布に従っていると考えることができます。試しに、少しだけ観測値を偏らせてP値がどうなるかをみてみます。
# 少しだけ偏った観測されたサンプル
observed = [450, 550, 450, 550, 450, 550]
result = stats.chisquare(observed, expected)
print(f'χ2={result.statistic} P-Value ={result.pvalue}')
直感的には、観測値がたった50ずれただけなので、このような観測値はありえなくはなさそうという気持ちになりますが、実際には0.00001の確率でしか発生しない奇跡であり、おそらく観測値は期待した分布とは違う分布から生成されたであろうであろうことがわかります。
ではこの適合度の$\chi^2$検定を使って、MCレイトレーシングのpdf関数とsample関数がちゃんと一致しているのかを確認してみましょう！
pdfとsampleを$\chi^2$検定する¶
pdf関数やsample関数の引数やサンプルは連続変数であり、カテゴリ変数ではないので、これを範囲を区切って、カテゴリ変数にします。つまりヒストグラムのようにすればよく、これはnp.histgram()を使うと変換ができます。
全てをまとめると次のような手順を踏みます。
- sample関数を用いて、いくつかのサンプルする。
- サンプルした結果をヒストグラムにする。
- そのヒストグラムの区間での理論確率を求めるためにpdf関数を呼ぶ。
- 2と3の結果を用いて$\chi^2$検定を行う。
# 定義域が0から1であることを仮定
def chiSquareTest(pdf, sample, n, bins):
u = np.random.random(n)
samples = sample(u)
observed, _ = np.histogram(samples, bins=bins, range=(0.0,1.0))
ibins = 1.0/bins
# TODO: bcenter = (brange[1:] + brange[:-1]) / 2 とかでできない？
x = np.linspace(ibins,1.0-ibins, bins)
expected = pdf(x) * n / bins
expected *= n/np.sum(expected)
#print(observed)
#print(expected)
result = stats.chisquare(observed, expected)
return result.pvalue
では完成した適合度の$\chi^2$検定を使って、sin関数のpdf関数と、偽物のpdf関数で結果がどのように変わるかをみてみましょう！
# 本物のpdfのP値
print(f'P-Value(sin_pdf): {chiSquareTest(sin_pdf, sin_sample, 10000, 50):10f}')
# 偽物のpdfのP値
print(f'P-Value(sin_like_pdf): {chiSquareTest(sin_like_pdf, sin_sample, 10000, 50):10f}')
P値が全く違うものになりました！つまり二次関数はsin関数ではなく偽物だということが明らかになりました。サンプル数とビンの数は上のプロットに使ったものと同じであることに注目してください。つまり目で見ても違いがあまり明確ではない場合でも、$\chi^2$検定を用いれば、明確な差として現れるということになります。驚きですね。
まとめ¶
適合度の$\chi^2$検定を用いて、MCレイトレーシングのpdf関数とsample関数をvalidateする方法について述べました。
ぜひ自作のMCレイトレーサーに$\chi^2$検定を組み込んで、コードの堅牢性を高めましょう！
# cf. https://www.anderswallin.net/2016/05/scipy-stats-chi2-ppf-without-scipy/
# 第1種不完全ガンマ関数
def lower_gamma(s,x):
g = 0
last_g = 1.0
done = False
tol = 1.0e-6
k=0
while not done:
top = pow(x, s)*math.exp(-x)*pow(x,k)
bot = np.prod( [float(s+j) for j in range(k+1) ] )
dg = float(top)/float(bot)
if dg == float("Inf"):
break
g += dg
k += 1
if k>100:
if g==0:
break
delta = abs(dg/g)
if delta == float("Inf"):
break
if delta < tol:
done = True
last_g = g
return g
# χ二乗分布のCDF
def chi2_cdf(x, k):
t = lower_gamma(k/2.0, x/2.0)
s = math.gamma(k/2.0)
return lower_gamma(k/2.0, x/2.0) / math.gamma(k/2.0)
# ユニバーサル関数化
chi2_cdf = np.frompyfunc(chi2_cdf,2,1)
# numpyを用いたchi2
def chisquare_numpy(observed, expected):
n = len(observed)
k = n - 1
n_p = np.array(expected)
res2 = (np.array(observed) - n_p) ** 2
chi2 = np.sum(res2/n_p)
ipv = chi2_cdf(chi2, k)
class Result:
statistic = chi2
pvalue = 1.0 - ipv
return Result()
# 観測される各カテゴリのサンプル数
observed = [468, 494, 532, 494, 517, 495]
# 理論上の各カテゴリのサンプル数。どのカテゴリも等しい数サンプルされることを期待する。
expected = [500, 500, 500, 500, 500, 500]
# scipy.statsを使って$\chi^2$検定をする
result_scipy = stats.chisquare(observed, expected)
result_mine = chisquare_numpy(observed, expected)
print(f'scipy χ2={result_scipy.statistic} P-Value ={result_scipy.pvalue}')
print(f'mine χ2={result_mine.statistic} P-Value ={result_mine.pvalue}')
# htmlに変換してダウンロード
import google.colab as colab
from google.colab import files
colab.drive.mount('/content/drive/')
%cd "/content/drive/MyDrive/colab/"
!jupyter nbconvert --to html "chi2.ipynb"
files.download('chi2.html')