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?

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