I. M. Sobol – The Monte Carlo Method – Section 08 – Calculation of the 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. 

One of widely used methods that we know about the calculation of the definite integral is the method of infinitesimal rectangles. We know that we get a more and more precise answer as the width of the rectangles decreases. In this page, we will look at the Monte Carlo Method of calculation of definite integrals. The Algorithm is outlined below.

  • 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 ξ
  • Use this ξ tow calculate two quantities
    • PDF[ ξ ]
    • [ ξ ]
  • Define η as G[ξ]PDF[ξ]
  • Generate several of these η1, η2, η3, η4, …. ηn
    • The distribution of the mean of these 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 η 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]

End of the Post 😉


.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.