Setup

We consider the function \(\varphi(x, y ) = 0.5+2sin\Big(cos(x)\sqrt{1+x^2+y^2}\Big)\) on the support \(S= [-1, 0]\times[-3, -2]\).

3D plot of function \(\varphi(x,y)\) and its contourplot


Notice that, as required, this function

We use 3 approaches to calculate the integral \(\int_{-3}^{-2}\int_{-1}^{0} \varphi(x,y) dxdy\).

Standard numeric calculation (using the Gaussian quadrature) yields the value of integral \(1.85502\).

1. Naïve Monte Carlo

We simulate samples of \(n=10^3\) i.i.d. points: \(X_i\sim U_{[-1, 0]}, Y_i\sim U_{[-3, 2]}, W_i\sim U_{[0, K]},\) where \(K = \max\limits_{(x, y)\in S}\varphi(x, y) = 2.5\).

We estimate integral as \(I_1 = K\cdot Vol(S)\frac{\sum_i^n(\varphi(X_i,Y_i)>W_i)}{n}\)

I1 <- (b-a)*(d-c)*K*sum(f(U,V) > W)/n

We plot the histogram of 100 integral values calculated using this method. The values are concentrated around the mean = \(1.71174\), which is quite lower than the target value.

2. Basic Monte Carlo

We simulate samples of \(n=10^3\) i.i.d. points: \(X_i\sim X_{[-1, 0]}, Y_i\sim U_{[-3, 2]}\) We estimate integral as \(I_2 = Vol(S)\cdot\frac{\sum_i^n(\varphi(X_i,Y_i))}{n}\)

I2 <- (b-a)*(d-c)*mean(f(U,V))

Again, we plot the histogram of 100 integral values calculated using this method. The values are concentrated around the mean = \(1.856057\), which is quite better than the result of previous method.

3. Importance Sampling

a) Independent margins

Now, we suppose that \(X_i, Y_i, i\in \{1,...,n\}\) are independent but not uniform on \([-1, 0 ]\) and \([-3, -2]\), respectively.

We suggest \(f(x) = -x+\frac{1}{2}\) to be the density of \(X\) and \(g(y)=y+3.5\) to be the density of \(Y\)

We make sure \(\int_{-1}^{0}f(x)dx=1\), \(\int_{-3}^{-2}g(y)dy=1\).

Next, the couple \((X, Y)\) has joint distribution \(p(x,y) = f(x)\cdot g(y)=(-x+\frac{1}{2})\cdot(y+3.5)\) and make sure \(\iint_S p(x,y) dS=1\)

3d plots of functions \(\varphi(x,y)\) (in rainbow) and the joint distribution p(x,y) (in green)

Remark: Notice that \(p(x, y)\) (in green) is a plot of joint density, hence it must integrate to 1, and it has different scale than \(\varphi(x,y)\). It is important that its slope is similar to that of \(\varphi(x,y)\) but its scale may differ.

We simulate sample of \(n=10^3\) couples \((X, Y)\sim p\) using inversion method.

To use inversion method, we calculated c.d.f. and their inverses, respectively:

\(F(x)=0.5(-x^2+x+2)\), \(F^{-1}(u)=0.5(1-\sqrt{9-8u})\) \(G(y) = 0.5(y^2+7y+12)\), \(G^{-1}(v)=-3.5+0.5\sqrt{8v+1}\)

We estimate integral as \(I_{3a}=\frac{1}{n}\sum_i^n\frac{\varphi(X_i,Y_i)}{p(X_i,Y_i)}\)

U <- runif(n,0,1)
V <- runif(n,0,1)
X <- Finv_u(U)
Y <- Ginv_v(V)
I3 <- 1/n*sum(f(X,Y)/p(X,Y))

The histogram of 100 integral values shows they are concentrated around the mean = \(1.85554\), which is closer to the target value and better than the result of previous method.

b) Dependent margins

Finally, we add dependence to the distribution of \((X,Y)\) using Frank copula. We used Frank copula with \(\theta = 0.3\) to make sure that couples where \(X\) and \(Y\) are both large or both low would occur more frequently.

We preserve the marginal distributions of \(X\) and \(Y\) to be the same as in the previous paragraph, \(f(x) = -x+\frac{1}{2}, g(y)=y+3.5\)

We set the joint distribution \(p_{Fr}(x,y)=f(x)g(y)c(F(x),G(y))\), where \(c(u,v) = \frac{\partial^2 C(u,v)}{\partial u\partial v}= \frac{e^{-\theta u}(-\theta e^{-\theta v)}(e^{-\theta}-1)}{(e^{-\theta}-1+(e^{-\theta u}-1)(e^{-\theta v}-1))^2}\)

The function \(\varphi(x,y)\) and the joint dencity \(p_{Fr}(x,y)\) look very similar. The closet the ration between them is to constant, the better importance sampling works. Let’s check it.

We generate \(n=10^3\) couples \((X,Y)\sim p_{Fr}\) as follows.

We generate n i.i.d. couples \((U_i,V_i)\sim C\), and we set \((X,Y)=(F^{-1}(U), G^{-1}(V))\).

In order to simulate \((U_i,V_i)\sim C\), we calculated the inverse function of \(F_u(z)=\frac{\partial C(u,v)}{\partial u}\): \(F_u^{-1}(z)=-\frac{1}{\theta} log\Big(1+\frac{e^{\theta u}z(e^{-\theta}-1)}{1-(1-e^{\theta u})z}\Big)\)

U <- runif(n) 
Z <- runif(n)
V <- dC_frankdu_inv(U,Z)
X <- Finv_u(U)
Y <- Ginv_v(V)
I3b <- 1/n*sum(f(X,Y)/p(X,Y))

Comparison of convergence rate

Finally, we compare the convergence of discussed methods as \(n\to \infty\). We plot the values of integral calculated using different methods as functions of \(n\).