# Monte Carlo Integration

I need to find the center of mass of a bar of length 12 using basic Monte Carlo integration. The mass distribution $$m(x) = .06x^2-3x+36$$ and the definition of the center of mass is

$$x_{com} = \frac{\int_0^{12}xm(x)dx}{\int_0^{12}m(x)dx}$$

Therefore, in the Wolfram Language(non Monte Carlo) we have

m[x_]:=.06x^2-3x+36
Integrate[x m[x],{x,0,12}]/Integrate[m[x],{x,0,12}]
= 4.68966


The Monte Carlo Method employs the fact that $$\lim_{n\rightarrow\infty}\frac{1}{n}\left(\sum_{i=1}^nf(x_i)\right) =$$ The average of the function on the domain $$[a,b]$$ where $$x_i$$ is a uniformly distributed random real such that $$a.

$$f_{avg} = \frac{1}{b-a}\int_a^bf(x)dx$$ and so we can draw the following conclusion:

$$\int_a^bf(x)dx\approx\frac{b-a}{n}\left(\sum_{i=1}^nf(x_i)\right)$$

Therefore, I write

n=100000;
a=0;
b=12;
xi:=RandomReal[{a,b}]
mcNumerator = (1/n)Sum[xi m[xi],{i,1,n}];
mcDenominator = (1/n)Sum[m[xi],{i,1,n}];
MC = mcNUmerator/mcDenominator


The output of this little program converges to something near 6 for large $$n$$. What am I missing here?

Edit I see a potential problem in the mcNumerator where each of those $$x_i$$ are different.

• xi generates a new random number everytime it appears. So xi m[xi] actually generates two random numbers. But Monte-Carlo integrationswants to have the same random number both times. Also, you need the factor 1/n in both the numerator and denominator, no? Feb 2 at 18:20
• I saw the multivaluedness in mcNumerator after I posted and 1/n was a typo. Fixed. Thanks. I will rewrite for get xi the same in xi m[xi] and test.
– JEM
Feb 2 at 19:01
• Just in case, surely you are aware of the Method -> "MonteCarlo" option on NIntegrate? Feb 3 at 11:39

The problem in your original code is that it generates different random values whenever you access the symbol xi (which is a function without arguments), so there are two different xi values in the xi * m[xi] - thus you actually computed two separated integrals for the numerator; more precisely, instead of computing

Integrate[x*m[x],{x,a,b}]


you computed

Integrate[x,{x,a,b}] * Integrate[m[x],{x,a,b}] / (b-a)


The numerator/denominator was then

Integrate[x,{x,0,12}] / 12


which is indeed 6.

The solution you proposed may be slightly inefficient - by appending one by one thousands (if not millions) of elements to the list there might be several reallocations and moving of data.

Still, your original version can be corrected by minor alterations:

m[x_]:=.06x^2-3x+36;
n=10000;
a=0;
b=12;
xi=RandomReal[{a,b}, n];
mcNumerator = ((b-a)/n)Total[xi* m[xi]];
mcDenominator = ((b-a)/n)Total[m[xi]];
MC = mcNumerator/mcDenominator


Of course, you may drop the (b-a)/n in front the two variables - I used them to check if each integral is computed correctly.

• Very nice answer. Your explanation of why the original code converged to 6 was very helpful.
– JEM
Feb 3 at 19:50

For something different:

In the following 2 rectangles that contain the regions of interest in numerator and denominator are randomly populated and fraction ‘under the respective curves’ is calculated and scaled:

m[x_]:= 0.06 x^ 2 - 3 x + 36
r1 = ReIm@ RandomComplex [ {0, 12 + 130 I},100000];
r2 = ReIm@RandomComplex [ {0, 12 + 36 I},100000];
al = (True / . CountsBy[r1, 0<#[]<[]m[#[1] &]) 12 × 130 / 100000.
a2 = (True / . CountsBy[r2, 0<#[] <m[#[ 1] &]) 12 × 36 / 100000.
a1 / a2


This yielded 4.687

Visualising:

ListPlot[GroupBy[r1, 0 < #[] <#[]m[#[]]&]]
ListPlot[GroupBy[r2, 0 < #[] <m[#[]] &]] I got this sorted out using a Do loop of all things:

m[x_] := .06 x^2 - 3 x + 36
data = {};
n=10000;
Do[{x = RandomReal[{0, 12}], AppendTo[data, {m[x], x m[x]}]},n]
t = Total[data]
t[]/t[]


This program converges to 4.68966

• Functional way to write this is something like With[{m = .06 #^2 - 3 # + 36 &}, Total[#2]/Total[#1] & @@ Transpose[{m[#], # m[#]} & /@ RandomReal[{0, 12}, 100000]]]. Maybe the use of anonymous functions (... &) a bit excessive on this, but using Do (especially in combination with AppendTo) is something most people working with Mathematica rather avoid, using Map, Table, Array and such instead. Feb 3 at 12:14
• ... or even With[{m = .06 #^2 - 3 # + 36 &}, Divide @@ Total[{# m[#], m[#]} & /@ RandomReal[{0, 12}, 100000]]]. Feb 3 at 12:21
• Or even Divide @@ Total[{#, 1} (.06 #^2 - 3 # + 36) & /@ RandomReal[{0, 12}, 100000]]. Okay, enough obfuscation is enough. :) Feb 3 at 12:34
• @ kirma I suspected the existence of an efficient functional approach to this but couldn't quite get there myself. Thank you. In your opinion, what is the best resource available that will take a traditional programmer and aid them in developing the kind of functional programming intuition you have displayed here?
– JEM
Feb 3 at 19:42
• It's an interesting question, I believe some of the old, and highest-voted answers here on Mathematica SE might dive more into it. I'm a bit biased on the subject since I learned to appreciate functional programming almost three decades ago on the introductory university CS course, and it stuck. I must admit that I was very much an imperative programming guy before that (but so was the hobbyist microcomputer world...). Feb 4 at 3:36
m[x_] := 0.06 x^2 - 3 x + 36
NIntegrate[x m[x], {x, 0, 12}]/NIntegrate[m[x], {x, 0, 12}]
n = 50000;
a = 0;
b = 12;
xi = Table[RandomVariate[UniformDistribution[{a, b}]], {i, 0, n}];
mcNumerator = Sum[xi[[i]] m[ xi[[i]] ], {i, 1, n}];
mcDenominator = Sum[m[xi[[i]]], {i, 1, n}];
MC = mcNumerator/mcDenominator


I apologize for not knowing how to enter Latex code to without triggering the "that's code; it must be indented" nagger.

m=Function[x, 0.06x^2-3 x+36]
x=RandomReal[12,10^4];
arm=Total[Map[# m[#]&,x]] ==>
mass=Total[Map[m,x]]
arm/mass ==> 4.70536
ClearAll[x]
\[Integral]x m[x]\[DifferentialD]x ==> 18 x^2-x^3+(3 x^4)/200
\[Integral] m[x] \[DifferentialD]x ==> 36 x-(3 x^2)/2+x^3/50
\*SubsuperscriptBox[$$\[Integral]$$, $$0$$, $$12$$]$$x\ m[x] \[DifferentialD]x$$\) ==> 1175.04
\*SubsuperscriptBox[$$\[Integral]$$, $$0$$, $$12$$]\ $$m[x] \[DifferentialD]x$$\) --> 250.56
1175.04/250.56 ==> 4.68966