2
$\begingroup$

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]);
            }
        }
    }
    answer = I/bin;
    cout << answer << endl;
}
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – A.S.
    Nov 13 '15 at 6:44
  • $\begingroup$ 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. $\endgroup$
    – A.S.
    Nov 13 '15 at 6:59
  • $\begingroup$ 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. $\endgroup$
    – Collaptic
    Nov 13 '15 at 7:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.