This is homework

(1)

已知： $$X_1, X_2, …, X_n \overset{\text{iid}}{\sim}p(x)$$

計算：

$$ E( \hat{I}_M)=E\left[\frac{1}{n} \sum^n_{i=1} \frac{f(X_i)}{p(X_i)} \right]=\frac{1}{n}E\left[ \sum^n_{i=1} \frac{f(X_i)}{p(X_i)} \right] $$

對於每個獨立的 $X_i$ ，我們只要計算： $$E\left[\frac{f(X_i)}{p(X_i)} \right]$$

因此： $$ E\left[\frac{f(X)}{p(X)} \right] = \int^b_a\frac{f(x)}{p(x)}p(x)dx =\int^b_af(x)dx = I $$

可知 $$E\left[\frac{f(X_i)}{p(X_i)} \right] =I,　\forall i $$

所以 $$ E(\hat{I}_M) =\frac{1}{n}\sum^n_{i=1}I=I $$

(2)

計算變異數 $$Var(\hat{I}_M)=E\left[(\hat{I}_M-I)^2\right]$$

因為 $$ \begin{aligned} Var(\widehat{I}_M) &= Var\left(\frac{1}{n} \sum_{i=1}^{n} \frac{f(X_i)}{p(X_i)}\right) = \frac{1}{n}Var\left(\frac{f(X)}{p(X)}\right) \\ &= \frac{1}{n}\left(E\left[\left(\frac{f(X)}{p(X)}\right)^2\right]-I^2\right) \end{aligned} $$

已知 $$E\left[\left(\frac{f(X)}{p(X)}\right)^2\right] < \infty$$

所以當 $n \to \infty$ 時 $$Var(\hat{I}_M) \to 0$$

因此 $$Var(\hat{I}_M)=E\left[(\hat{I}_M-I)^2\right]\to 0$$

即 $$\hat{I}_M \overset{L^2}{\to} 0$$

(3-a)

我們計算 $\hat{I}_M$的變異數，由(2)可知 $$ Var(\hat{I}_M)=\frac{1}{n}\left(E\left[\left(\frac{f(X)}{p(X)}\right)^2\right]-I^2\right) $$

其中 $$ \tag{1} E\left[\left(\frac{f(X)}{p(X)}\right)^2\right]=\int^1_0\frac{f(x)^2}{p(x)}dx $$

$$ \tag{2} I = E\left[\frac{f(X)}{p(X)} \right] = \int^b_a\frac{f(X)}{p(X)}p(x)dx $$

$$ \Rightarrow I= \int^1_0(x^{-1/3}+\frac{x}{10})dx \approx 1.55 $$

當 $p_1(x)=1$ 時，$x \in (0, 1)$，則 $$ \begin{aligned} E\left[\left(\frac{f(X)}{p_1(X)}\right)^2\right] &=\int^1_0f(x)^2dx \\ &=\int^1_0(x^{-1/3}+\frac{x}{10})^2dx \\ &\approx 3.1233 \end{aligned} $$

因此 $$ Var(\hat{I}_M)=\frac{1}{n}(3.1233-1.55^2)=\frac{0.7208}{n} $$

當 $p_2(x)=\frac{2}{3}x^{-1/3}$ 時，$x \in (0, 1]$，則 $$ \begin{aligned} E\left[\left(\frac{f(X)}{p_2(X)}\right)^2\right] &=\int^1_0\frac{f(x)^2}{p_2(x)}dx \\ &=\int^1_0\frac{(x^{-1/3}+\frac{x}{10})^2}{\frac{2}{3}x^{-1/3}}dx \\ &= 2.4045 \end{aligned} $$

因此 $$ Var(\hat{I}_M)=\frac{1}{n}(2.4045-1.55^2)=\frac{0.002}{n} $$ 由上述可知， $p_2(x)$ 會使變異數變小

import numpy as np

def f(x):
    return x**(-1/3) + x/10
def p1(n):
    return np.random.uniform(0, 1, n)
def p2(n):
    u = np.random.uniform(0, 1, n)
    return u**3
def monte_carlo_int(n, sample_f, p_f):
    X = sample_f(n)
    weights = f(X)/p_f(X)
    return np.mean(weights)

n_samples = 5000
n_test = 1000
estimate_p1 = np.zeros(n_test)
estimate_p2 = np.zeros(n_test)

for i in range(n_test):
    estimate_p1[i] = monte_carlo_int(n_samples, p1, lambda x:1)
    estimate_p2[i] = monte_carlo_int(n_samples, p2, lambda x: (2/3)*x**(-1/3))

var_p1 = np.var(estimate_p1, ddof=1)
var_p2 = np.var(estimate_p2, ddof=1)

print(f'The estimator variance by Monte-Carlo of p1(x) is: {var_p1: .6f}')
print(f'The estimator variance by Monte-Carlo of p2(x) is: {var_p2: .6f}')

The estimator variance by Monte-Carlo of p1(x) is:  0.000139
The estimator variance by Monte-Carlo of p2(x) is:  0.000000

(3-c)

為了讓 $g(x)$ 包含 $f(x)$，我們選擇一個常數 $\alpha$使得 $$e(x)=\alpha g(x) \ge f(x),　\forall x$$

而我們定義隨機變數 $N$ 為成功產生一個樣本 $X,　(X\in g(x))$ 所需的嘗試次數，意即

前 $N-1$ 次失敗，則 $U > f(x)/\alpha g(X)$
第 $N$ 次成功，則 $U \le f(x)/\alpha g(X)$

這表示 $$ P(N=k)=P(前k-1次失敗) *P(第k次成功)$$

因此，對於任意一次嘗試，成功的機率為 $$ p = \int^{\infty}_0g(x)\frac{f(x)}{\alpha g(x)}dx = \frac{1}{\alpha}\int^{\infty}_0f(x)dx =\frac{1}{\alpha} $$

由此可知每次嘗試成功的機率為 $\frac{1}{\alpha}$，而失敗的機率為 $1-p = 1-\frac{1}{\alpha}$

最後計算 $P(N=k)$，即前 $k-1$ 次失敗，第 $k$ 次成功的機率 $$P(N=k)=(1-p)^{k-1}p$$

代入 $\frac{1}{\alpha}$ $$P(N=k)=\left(1-\frac{1}{\alpha}\right)^{k-1}\frac{1}{\alpha}$$

可得 $$ N \sim Geo(p)$$

(3-d)

由幾何分布的 p.m.f. 可知 $$P(n=K) = (1-p)^{k-1}p,　k=1, 2, 3,…$$ 計算期望值 $$ \begin{aligned} E(N) &= \sum^\infty_{k=1}k (1-p)^{k-1}p = p\sum^\infty_{k=1}k (1-p)^{k-1} \\ &= p*\frac{1}{p^2}= \frac{1}{p} \end{aligned} $$

代入 $p=\frac{1}{\alpha}$ 可得 $$E(N)=\alpha$$

This is homework#

(1)#

(2)#

(3-a)#

(3-c)#

(3-d)#

This is homework

(1)

(2)

(3-a)

(3-c)

(3-d)