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?

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.
