Monte Carlo simulations for arbitrary functions

Question

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?

Answer 1

Draw $n$ pairs $(x,y)$, iid uniformly distributed in the unit square. Count how many of these pairs satisfy $y<x^2$, let this number be $k$. Then $$\mathbb P(Y<X^2) = \int_{[0,1]^2} \mathbb I_{y<x^2}\,\text d(x,y) = \int_0^1 x^2\,\text dx\approx\frac{k}{n}.$$

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).

Answer 2

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.

Answer 3

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

Addendum 2: About extending MC to multiple dimensions.

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.