# How to use Monte Carlo Integration for multi-dimensional integration

I'm having trouble implementing Monte Carlo for multi-dimensional integration in C++, and looking for some guidance.

Here is the function I'm suppose to integrate. $$I_n = \int_{0}^1 dx_1 \int_{0}^1 dx_2 ...\int_{0}^1 dx_n \sqrt (x_1\times x_2 \times...\times x_n)$$

I'm suppose to consider any integrand of n variables $f(\mathbf x)$, assuming each component $x_i$ of $\mathbf x$ is restricted to $x_i \in [0,1]$ and assuming $\;$ $f \in [0,1]$. I'm suppose to generate random samples for $x_i$ and create the resulting vector $\mathbf x^{'}$. I'm also tasked with generating random rejection samples for $a$. Finally, I am to determine if $a < f(\mathbf x^{'})$, and add 1 to the bin (i.e. under the curve). How do I do this with Monte Carlo Integration for n dimensions? By the way, m is the number of divisions along a given axis, $m^n$ is the number of evaluations of the integrand.

Here is my attempt at the problem, however I feel as if I'm going about this the wrong way. (I'll provide my code for my first 2 dimensions to keep it simple).

Edit: When I do n = 3, I get the answer = 0.422095

int main(int argc, char *argv[])
{
double size = pow(m,n);
vector<double> x1(size);
vector<double> x2(size);
vector<double> a(size);
double bin = 0;
double I=0;
for (int i = 0; i<size; ++i)
{
x1[i] = (double)rand()/(double)RAND_MAX;
x2[i] = (double)rand()/(double)RAND_MAX;

for (double i = 0; i<size; ++i)
{
a[i] = (double)rand()/(double)RAND_MAX;
}

if (n == 1)
{
for (double i = 0; i < size; ++i)
{
if (a[i] < sqrt(x1[i]))
{
bin = bin + 1;
I = I + sqrt(x1[i]);
}
}
}
else if (n == 2)
{
for (double i = 0; i < size; ++i)
{
if (a[i] < sqrt(x1[i]*x2[i]))
{
bin = bin + 1;
I = I + sqrt(x1[i]*x2[i]);
}
}
}
}

• 1. Generate $x=(x_1,\dots,x_n)=(U(0,1))^n$ 2. Generate $a=U(0,1)$. 3. If $f(x)<a$, increase count. 4. Repeat.
– A.S.
Nov 13 '15 at 6:44
• Also, for large $n$, CLT gives $f(x)\approx Log-Normal (-N,N)$, hence $I_n\approx exp(-\frac {3N}8)$ and you can use this to check your computational work.
– A.S.
Nov 13 '15 at 6:59
• Thanks for the help, for your equation exp(-3N/8), N is referring to the number of dimensions? I assume so, because I've played around with my code a bit more, and I believe I have it working now. When I do 3 dimensions, I get an answer of around .3056, which is close to 0.3246 when I use N=3 in your equation. Nov 13 '15 at 7:10