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?

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.
