# Monte Carlo simulations for arbitrary functions

I'm familiar with MC methods for approximating PDF integrals. But in this question, I'm curious how we might adapt these methods for other problems. For example evaluating $$\int_{0}^{1} x^2 dx$$ . I choose this function because evaluating the integral analytically is trivial; I'm simply curious how one would design an MC simulation to approximate the value that could be found analytically.

Edit @Periwinkle, I found the below code snippet (cleaned it up for ease of readability) and posted below.

def mc_int(upper, lower, size, func):
uniform_samples = np.random.uniform(low=lower, high=upper, size=size)
transformed_samples = func(uniform_samples)
expected_value = np.average(transformed_samples) * (upper - lower)
return expected_value

def f(x):
return x**2

mc_int(lower=0, upper=1, size=1000, func=f)
>>>
0.3333


Based on your comment, I don't grasp why scaling by upper - lower is necessary. Could you explain?

• In R, use function to define the specific function, then integrate (with your function as integrand and appropriate limits) to get the integral. Commented Oct 5, 2021 at 16:23
• On scaling by upper - lower: Try finding your integral for a range which is not of width $1$, for example $\int\limits_{0.4}^{0.41} x^2\,dx$. You will get a lot of random values between $0.4^2=0.16$ and $0.41^2=0.1681$ with an average of around $0.164$. But that is far too high for the integral itself which, considering bounding rectangles for the area under the curve, must be between $0.16(0.41-0.4)$ and $0.1681(0.41-0.4)$ because of the narrowness of the integral. Hence the need to scale by upper - lower. The correct answer is about $0.001640333$ Commented Oct 6, 2021 at 8:16
• @Henry, I don't know how I missed this! If $E(x)$ is the height, we still need a width to return an area. Thanks for the help Commented Oct 6, 2021 at 15:11

Draw $$n$$ pairs $$(x,y)$$, iid uniformly distributed in the unit square. Count how many of these pairs satisfy $$y, let this number be $$k$$. Then $$\mathbb P(Y

R code:

xx <- seq(0,1,by=.01)
plot(xx,xx^2,type="l",lwd=3)

n_sims <- 1e4
set.seed(1)
sims <- cbind(runif(n_sims),runif(n_sims))
index <- sims[,2]<sims[,1]^2
points(sims,pch=19,cex=0.4,col=c("red","black")[index+1])

sum(index)/length(index)


This is a variation of a well-known exercise about approximating $$\pi$$ (actually $$\frac{\pi}{4}$$, by using a quarter-circle in the unit square).

• Would this work on say multivariable functions, say $z = f(x,y)$? My gut instinct says yes, but that the number of points to sample would grow exponentially with variables in the system; seems like the curse of dimensionality would be applicable here. Commented Oct 5, 2021 at 16:46
• This was my thought, but how would you know that a function is contained within the unit square? This seems particularly problematic for funky distributions in many dimensions.
– Dave
Commented Oct 5, 2021 at 17:15
• relevant wiki page: en.wikipedia.org/wiki/Monte_Carlo_integration
– eps
Commented Oct 6, 2021 at 4:30
• @jbuddy_13: yes to both. It would work for $N$-variate functions (you would need to draw $(N+1)$-tuples, the first $N$ in the domain of the function, the last one between the minimum and the maximum of $f$). And as you write, the number $n$ of samples required for a given quality of the approximation grows very fast with the dimensions. But then again, the same holds for many other numerical approximation methods. Commented Oct 6, 2021 at 6:04
• @Dave: yes, you are completely right. Such an MC integration depends crucially on (1) knowing the minimum and the maximum of the function over the integration range (kind of hard for functions that are unbounded but integrable), and (2) being able to quickly evaluate the function so we know whether a sample is "below" or "above" the function. Commented Oct 6, 2021 at 6:06

The representation of$$\mathfrak I = \int_0^1 x^2\,\text dx$$as an expectation of a random variable is quite open, in that the choice of a Uniform (0,1) variable $$U$$ such that$$\mathfrak I = \mathbb E[U^2]$$is not the only choice. For any strictly positive probability density $$f$$ over $$(0,1)$$, the representation$$\mathfrak I = \int_0^1 x^2\,f^{-1}(x)\,f(x)\,\text dx$$ holds, meaning that $$\mathfrak I = \mathbb E^f[X^2\,f^{-1}(X)]$$can be exploited in a Monte Carlo scheme: $$\mathfrak I = \lim_{n\to\infty}\frac{1}{n}\sum_{i=1}^n X_i^2\,f^{-1}(X_i)\qquad X_i\sim f(x)$$The (formally) optimal choice of $$f$$ is the Beta $$\mathcal B(2,1)$$ density, in that the resulting Monte Carlo approximation has variance zero.

• Of course the formally optimal choice requires you to already know the answer in order to set up the implementation.
– Ian
Commented Oct 6, 2021 at 13:22
• Of course (and this is why I used the term formally!, with an unintentional typo) Commented Oct 6, 2021 at 13:36

My purpose here is to show a Riemann approximation for $$\int_0^1 x^2\, dx = 1/3$$ with enough rectangles to approximate the integral. Then to do a Monte Carlo integration in which uniformly chosen points in the interval of integration are substituted for centers of bases of rectangles.

# Riemann approx with m rectangles
m = 1000; a = 0; b = 1
w = (b-a)/m  # rectangle widths
d = seq(a+w/2,b-w/2, len=m) # centers
h = d^2  # rectangle heights
sum(w*h) # rectamg;e areas
[1] 0.3333332

# MC emulates Riemann
# with random uniform grid
m = 10^6;  a=0;  b=1
w = (b-a)/m   # "average width"
d = runif(m, a, b)
h = d^2
sum(w*h)
[1] 0.3332943


Addendum 1: per question in comment by @Dave: MC approximation of the density function of $$\mathsf{Beta}(2,2)$$ to verify it integrates to unity.

# Approx integration: BETA(2,2) density
m = 10^6;  a=0;  b=1
w = (b-a)/m
d = runif(m, a, b)
h = dbeta(d, 2, 2)  # BETA(2,2) PDF
sum(w*h)
[1] 1.000299  # aprx 1


Suppose we want to verify the probability $$0.3413^2 = 0.1165$$ in the unit square under a standard bivariate normal distribution. We put points uniformly at random in the unit square and sum their corresponding densities:

set.seed(1234)
m = 10^4; u1=runif(m); u2 = runif(m)
h = dnorm(u1)*dnorm(u2)
mean(h)                 # mc aprx
[1] 0.1163528
diff(pnorm(c(0,1)))^2   # exact
[1] 0.1165162


Also, we find the probability $$0.0677$$ of the standard bivariate normal distribution within the triangle with vertices $$(0,0), (0,1), (1.0).$$ We put points uniformly at random in the triangle, sum the corresponding values of the density, and multiply by the area $$1/2$$ of the triangle. [The exact value can be obtained by symmetry, using a 45-degree rotation.]

h.acc = h[u1+u2<=1]
.5*mean(h.acc)                # mc aprx
[1] 0.06768097
diff(pnorm(c(sqrt(1/2),0)))^2 # exact
[1] 0.06773003


Perhaps see this Q&A for additional MC integration methods.

• Does this rely on us knowing that the function is bounded by $0$ and $1$ on the $y$-axis?
– Dave
Commented Oct 5, 2021 at 17:18
• As presented, positive and bounded, Yes. But not necessarily by 1. See addendum with density of $\mathsf{Beta}(2,2),$ for which max is 1.5. Commented Oct 5, 2021 at 17:26