Throw a fistful of darts, blindfolded, at a square dartboard with a quarter-circle painted on it. Each dart lands somewhere at random. Count what fraction land inside the arc — and that fraction, multiplied by four, is π. No trigonometry, no infinite series, no clever formula: just counting hits. The same trick estimates any area or integral, including ones with no closed form at all — you fence off a region, rain random points into its bounding box, and the hit-rate times the box area is the answer. Watch a deterministic constant assemble itself out of pure chance, and watch the error shrink like 1/√N.
Put a quarter-circle of radius 1 inside the unit square \([0,1]\times[0,1]\). The square has area 1; the quarter-circle has area \(\pi/4\). Now scatter points uniformly at random across the square. Because the points are uniform, the probability that any one lands inside the quarter-circle equals the region's fraction of the total area — exactly \(\pi/4\). So if you throw \(N\) darts and \(N_{\text{in}}\) land inside (i.e. satisfy \(x^2+y^2\le 1\)), the ratio \(N_{\text{in}}/N\) estimates \(\pi/4\), and \(\pi \approx 4\,N_{\text{in}}/N\). Nothing about circles or angles is used — only counting and the definition of area as "how much of the whole a shape occupies."
The same idea integrates anything. To find \(\int_a^b f(x)\,dx\) for a non-negative \(f\), enclose the graph in a bounding box of known area \(A_{\text{box}}\), rain darts into the box, and count the fraction that land under the curve (the test is simply \(y \le f(x)\)). Then the area under \(f\) is \(A_{\text{box}}\cdot(N_{\text{in}}/N)\). This is why Monte Carlo shrugs at functions with no closed-form antiderivative, like \(e^{-x^2}\): it never needs a formula, only the ability to ask "is this point inside?". Swap the circle for a heart, a fractal, or a tangle of inequalities in a hundred variables and the recipe is unchanged.
What does it cost? Each dart is a Bernoulli trial with success probability \(p\) (the true area fraction). The count \(N_{\text{in}}\) is binomial, so the estimate \(\hat p = N_{\text{in}}/N\) has standard deviation \(\sqrt{p(1-p)/N}\). That is the crucial law: the typical error decays like \(\boldsymbol{1/\sqrt{N}}\). To halve the error you need four times as many darts; to win one more decimal digit you need a hundred times as many. On a log–log plot of error against \(N\) the points ride a line of slope \(-\tfrac12\) — you can watch that line assemble itself in mode ③. This slow convergence is Monte Carlo's price of admission, and it is why estimating π this way is charming but wildly inefficient compared to a series.
So why is the method everywhere — pricing options, simulating reactor neutrons, rendering film lighting, sampling Bayesian posteriors? Because that \(1/\sqrt{N}\) rate is blind to dimension. A grid method that places points every \(h\) apart needs \(h^{-d}\) of them in \(d\) dimensions — the "curse of dimensionality" makes it hopeless past a handful of variables. Monte Carlo's error depends on the variance of the integrand and the number of samples, not on \(d\). In ten or ten-thousand dimensions the darts still work, just as they do on our little dartboard. Randomness, it turns out, is not the enemy of precision — with enough samples it is a perfectly good way to compute the most deterministic quantities we know.