# 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? Commented Feb 2, 2023 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
Commented Feb 2, 2023 at 19:01
• Just in case, surely you are aware of the Method -> "MonteCarlo" option on NIntegrate? Commented Feb 3, 2023 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
Commented Feb 3, 2023 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];
a1 = (True /. CountsBy[r1, 0 < #[[2]] < #[[1]] m[#[[1]]] &])  12*130/
100000.;
a2 = (True /. CountsBy[r2, 0 < #[[2]] < m[#[[1]]] &])  12*36/100000.;
a1/a2


This yielded 4.687

Visualising:

ListPlot[GroupBy[r1, 0<#[[2]]< #[[1]] m[#[[1]]] &]]
ListPlot[GroupBy[r2, 0<#[[2]]<m[#[[1]]] &]]


Correction 2024-04-06

@Mariusz Iwanuik correctly identified that my original code did not run. My original post contained multiple typographical errors (related to OCR errors pasting from image of code on cloud). I have corrected and confirmed as screenshot demonstrates.

• The code dosen't work ? Commented Apr 6 at 9:14
• @MariuszIwaniuk thank you very much for identifying this. I hope I have corrected. I have stopped using OCR pastes because of the typographical errors. Thank you again. I have added a corrigendum. Commented Apr 6 at 9:51

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[[2]]/t[[1]]


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. Commented Feb 3, 2023 at 12:14
• ... or even With[{m = .06 #^2 - 3 # + 36 &}, Divide @@ Total[{# m[#], m[#]} & /@ RandomReal[{0, 12}, 100000]]]. Commented Feb 3, 2023 at 12:21
• Or even Divide @@ Total[{#, 1} (.06 #^2 - 3 # + 36) & /@ RandomReal[{0, 12}, 100000]]. Okay, enough obfuscation is enough. :) Commented Feb 3, 2023 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
Commented Feb 3, 2023 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...). Commented Feb 4, 2023 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