It should say f(xi)
f() is the function we are trying to integrate via the numerical monte carlo method, which estimates an integral (and its error) by evaluating randomly choosen points from the integration region.
Ref.
Your goal is to compute the integral of f from x1 to x2. For example, you may wish to compute the integral of sin(x) from 0 to pi.
Using Monte Carlo integration, you can approximate this by sampling random points in the interval [x1,x2] and evaluating f at those points. Perhaps you'd like to call this MonteCarloIntegrate( f, x1, x2 ).
So no, MonteCarloIntegrate does not "feed back" into itself. It calls a function f, the function you are trying to numerically integrate, e.g. sin.
Replace f(x_r) by f(x_r_i) (read: f evaluated at x sub r sub i). The r_i are chosen uniformly at random from the interval [x_1, x_2].
The point is this: the area under f on [x_1, x_2] is equal to (x_2 - x_1) times the average of f on the interval [x_1, x_2]. That is
A = (x_2 - x_1) * [(1 / (x_2 - x_1)) * int_{x_1}^{x_2} f(x)\, dx]
The portion in square brackets is the average of f on [x_1, x_2] which we will denote avg(f). How can we estimate the average of f? By sampling it at N random points and taking the average value of f evaluated at those random points. To wit:
avg(f) ~ (1 / N) * sum_{i=1}^{N} f(x_r_i)
where x_r_1, x_r_2, ..., x_r_N are points chosen uniformly at random from [x_1, x_2].
Then
A = (x_2 - x_1) * avg(f) ~ (x_2 - x_1) * (1 / N) * sum_{i=1}^{N} f(x_r_i).
Here is another way to think about this equation: the area under f on the interval [x_1, x_2] is the same as the area of a rectangle with length (x_2 - x_1) and height equal to the average height of f. The average height of f is approximately
(1 / N) * sum_{i=1}^{N} f(x_r_i)
which is value that we produced previously.
Whether it's xi or xr is irrelevant - it's the random number that we're feeding into function f().
I'm more likely to write the function (aside from formatting) as follows:
(x2-x1) * sum(f(xi))/N
That way, we can see that we're taking the average of N samples of f(x) to get an average height of the function, then multiplying by the width (x2-x1).
Because, after all, integration is just calculating area under the curve. (Nice pictures at http://hyperphysics.phy-astr.gsu.edu/Hbase/integ.html#c4.
x_r is a random value from the integral's range.
Substituting Random(x_1, x_2) for x_r would give an equivalent equation.