Simulation | Monte Carlo Method | Definite Integral

Note: In this section of the text, the author explains several related topics before getting the question. For this reason, I have chosen to not go with the Author's flow. 

• Choose the function for which the definite integral needs to be found. Let us call this \(G\)
  • Decide on the range for integration.
  • Let's call the limits as [a,b]
• Choose an arbitrary probability distribution function PDF that is defined in this interval [a, b]. This means that the following are implied
  • A valid CDF exists in the same range [a, b]
  • A valid Inverse CDF also exists in the same range [a, b]
  • what i mean to say is that the PDF has to be only confined to the range [a, b]
• Define a random variate off the PDF defined earlier
  • this means we would be using the Inverse CDF
  • Use this random variate to generate random number called \(\xi\)
• Use this \(\xi\) tow calculate two quantities
  • \(PDF(\xi\))
  • \(G(\xi)\)
• Define \(\eta = \frac{G(\xi)}{PDF(\xi)}\)
• Generate several of these \(\eta_{1}, \eta_{2}, \eta_{3}, \eta_{4} ... \eta_{n}\)
• The distribution of the mean of the above variables would be centered around the actual value of the integral

Simulation

• The function to be integrated is \(Sin(0.5 \pi x)\) in the range \([0, 2]\)
• The probability function chosen for this range is Triangular with range [0, 2] and the apex at 1
• The following plot shows the PDF, the CDF and the function \(Sin(0.5 \pi x)\)
• Note the following
  • the PDF and the function of the interest are have the same valid range
  • as in they are both defined from 0 to 2
• Since the triangular function is also a PDF, it is implied that the area under it is equal to one

After you have these functions setup, there are only few steps to the calculation of the integral (OR) more aptly the estimation of the integral. Using \(\eta\) as per the algorithm outlined above, we estimate the integral. The following series of images shows the integral calculated using various sample sizes. Observe how the estimation of the integral more and more precise when the sample size increases. The accuracy has not changed much but the precision greatly becomes better 😀

Code

    ClearAll[toIntegrate, probabilityFunction, cdfProbabilityFunction];

    toIntegrate[x_] := Sin[0.5 \[Pi] x]
    probabilityFunction[x_] :=
        Module[{constant = NIntegrate[-(r - 1)^2 + 1, {r, 0, 2}]}, Divide[-(x - 1)^2 + 1, constant]]
    probabilityFunction[x_] := PDF[TriangularDistribution[{0, 2}], x]
    cdfProbabilityFunction[x_] := Integrate[probabilityFunction[r], {r, 0, x}]
    inverseCDFProbabilityFunction[x_] := InverseCDF[TriangularDistribution[{0, 2}], x]

    Plot[
        {toIntegrate[x], probabilityFunction[x], cdfProbabilityFunction[x]}, {x, 0, 2}
        , PlotLegends -> "Expressions"
        , AspectRatio -> 0.5
        , ImageSize -> 788
        ,   Epilog -> {
                {Opacity[0.5], Blue, Dotted, InfiniteLine[{{1, 0}, {1, 1}}]},
                {Opacity[0.5], Dotted, Blue, InfiniteLine[{{0, 0.5}, {2, 0.5}}]}
            }
    ]

    Histogram[
      Mean[With[{\[Xi] = RandomVariate[TriangularDistribution[{0, 2}]]},
        toIntegrate[\[Xi]] / probabilityFunction[\[Xi]]] & /@
          Range@10] & /@ Range[10000], 200]

Reference Chapter from I.M. Sobol - The Monte Carlo Method

This topic is discussed on page 61 of the PDF below. Feel free to download the document for future references