Multiple Importance Sampling

Estimator:

\[F=\sum_i^n {1 \over n_i} \sum_j^{n_i} \omega_i (X_{i,j}) \frac{f(X_{i,j})}{p_i(X_{i,j})}\]

where

• $\sum_i^n \omega_i(x)=1$ whenever $f(x)\ne 0$
• $\omega_i(x)=0$ whenever $p(x)=0$

Expected value of estimator:

\[E[F]=\int_\Omega \sum_i^n {1\over n_i} \sum_j^{n_i} \omega_i(x) {f(x) \over p_i(x)} p_i(x) dμ(x)\]

注意这里将采样 $X_{i,j}$ 换成了随机变量x，这使得

\[\sum_j^{n_i} ω_i(x) f(x) = n_i \ ω_i(x) f(x)\]

因此上式可化简为：

\[E[F] = \int_\Omega \sum_i^n {1 \over n_i} n_i \ ω_i(x) f(x) dμ(x) = \int_\Omega \sum_i^n ω_i(x) f(x) dμ(x) = \int_\Omega f(x) dμ(x)\]

因此F可用于估算 $\int_\Omega f(x) dμ(x)$。注意 $F$ 估算的是整个积分，不是被积函数 $f(x)$。

Example

Constant weights

\[F={1 \over n_1} \cdot \omega_1 {f(X_{1,1})\over p_1(X_{1,1})} + {1 \over n_2} \cdot ω_2 {f(X_{2,1})\over p_2(X_{2,1})} + {1 \over n_3} \cdot ω_3 {f(X_{3,1}) \over p_3(X_{3,1})}\]

权重是常量： $\omega_i(x)=c_i$，

且满足：$\omega_1 + \omega_2 + \omega_3 = 1$，且 $n_i$ = 1。

这种情况，$F$ 是3个采样的加权平均。

The balance heuristic

\[\hat{\omega}_i(x)={n_i p_i(x) \over \sum_k^n n_k p_k(x)}\]

代入 $F$：

\[\hat{F} =\sum_i^n {1 \over n_i} \sum_j^{n_i} {n_i p_i(X_{i,j}) \over \sum_k^n n_k p_k(x)} {f(X_{i,j}) \over p_i(X_{i,j})} =\sum_i^n \sum_j^{n_i} {f(X_{i,j}) \over \sum_k^n n_k p_k(X_{i,j})} ={1 \over \sum_k^n n_k} \sum_i^n \sum_j^{n_i} {f(X_{i,j}) \over {\sum_k^n n_k p_k(X_{i,j}) \over \sum_k^n n_k }}\]

记：

\[N=\sum_k^n n_k\] \[c_k = {n_k \over \sum_k^n n_k} ={n_k \over N}\]

则：

\[\hat{F}={1 \over N} \sum_i^n \sum_j^{n_i} {f(X_{i,j}) \over \sum_k^n c_k p_k(X_{i,j})}\]

若记:

\[\hat{p}=\sum_k^n c_k p_k(X_{i,j}) = {\sum_k^n n_k p_k(X_{i,j}) \over \sum_k^n n_k}\]

则上式即表示为标准 Monte Carlo Estimator 的形式：

\[\hat{F}={1\over N} \sum_i^n \sum_j^{n_i} {f(X_{i,j}) \over \hat{p}(X_{i,j})}\]

上式也是最终用于计算的公式。仍然进行N次采样，但每次采样的pdf，需要根据对多个pdf进行加权平均。

伪代码

# Calculate N
for k in range(n):
    N += n_k

F = 0.0
for i in range(n):
    for j in range(n_i):
        X = Sample(i, j)
        # Calculate p
        p = 0.0
        for k in range(n):
            p += (n_k / N) * p_k(X)
        # Calc function
        f = func(X)
        # Accumulate sample
        F += f / p
F /= N