What does it mean when we say we have computed a number to a certain accuracy using a probabilistic algorithm?

Question

My intention is to ask a general question about probabilistic (Monte Carlo) algorithms. But to keep things simple, I will focus on a few specific examples.

Let me start the discussion with deterministic algorithms. Suppose I use some iterative algorithm to compute say the $\sqrt{2}$ (for example, using Newton's method) and I claim that if I run $n$ iterations of the algorithm, then the answer I obtain has an error of at most $\epsilon_n$. The precise meaning of this statement is that $$ |x_n - \sqrt{2}| \leq \epsilon_n,$$ where $x_n$ is the answer I get in the $n^{th}$ iteration.

Now suppose I have some iterative Monte Carlo algorithm. To take a simple example, suppose I want to compute the value of $\pi$ using the well known Monte Carlo algorithm: consider a square of side length $1$ and at each iteration, generate a pair of random numbers between $0$ and $1$. Check whether the point lies inside the circle quadrant or not. After $n$ iterations, calculate the fraction of points that are inside that circle quadrant. Call this number $\frac{\pi_n}{4}$.

$\textbf{Question 1:}$ What is a meaningful question one can ask about $\pi_n$? Saying something like $$ |\pi - \pi_n| \leq \epsilon_n $$
doesn't make any sense, because there is a small chance that $\epsilon_n$ could we very big (even if $n$ is very large).

$\textbf{Note:}$ To keep things simple let us assume that the random numbers we produce are "truly random" (whatever that means).

$\textbf{Question 2}:$ I now have a more specific question. It seems that using a Monte Carlo algorithm (called the Pivot Method), one can numerically compute the connective constant of a lattice (for this discussion it doesn't matter too much what connective constant means except that its a real number, so I won't bother defining it). Now consider the following statement:

The connective constant ($\mu$) of the square lattice upto two decimal places is $2.63$.

What does this statement mean? If that $2.63$ was obtained by an ordinary algorithm, then it would have meant $$|\mu-2.63| \leq 0.005. $$ But this $2.63$ was obtained by a Monte Carlo Algorithm. So its not clear to me what is actually meant by saying a number has been computed to some accuracy.

A discussion on connective constant (and their known numerical values) is available in the wikipedia page

http://en.wikipedia.org/wiki/Connective_constant

Answer 1

For Question 1, it seems you are looking for the notion of a confidence interval.

Answer 2

There is no universal standard for what an error estimate on a Monte Carlo algorithm means.

What you can say in general is this: when using a Monte Carlo algorithm to compute the value of some parameter $x$ whose true value is $x_\text{true}$, the algorithm gives a result $x_\text{MC}$ drawn from a probability distribution $\mathcal{P}_\text{true}$ parametrized by $x_\text{true}$. Given a specific result $x_\text{MC}$ and the distribution $\mathcal{P}_\text{true}$, one can determine the range of values for $x_\text{true}$ for which the likelihood ratio $\frac{\mathcal{L}_\text{true}}{\mathcal{L}_\text{MC}}$ exceeds a certain threshold. This range is the confidence interval that Bjørn mentioned. The Monte Carlo error estimate is a measure of the size of the confidence interval. (Incidentally, physical measurements with uncertainties are interpreted in the same way, as confidence intervals.)

Note that the size of the confidence interval depends on the chosen threshold for the likelihood ratio test, and in turn, so does the meaning of the error estimate. This is a choice made by the designer of the algorithm. When a result says $\mu=2.63$ to two decimal places, for example, it's expressing that the likelihood ratio exceeds a certain threshold over a range of size $\approx 0.01$ centered approximately on $2.63$, but that really tells you nothing unless you know what the threshold is.

In many common cases, the probability distribution is taken to be Gaussian (whether this is proven or just assumed for simplicity depends on the algorithm), and in that case the error estimate will usually be taken as the standard deviation of the distribution.

Answer 3

It means pretty much the same thing as when we say that we have measured a physical property to a certain accuracy. That is, in loose terms, it means that we are reasonably certain that the correct value lies close to the estimated value — and, if we're to be proper, we should also state how close, how certain we are of this, and how this certainty was estimated.

As Bjørn Kjos-Hanssen already noted, this (un)certainty is typically expressed by a confidence interval, which basically says that we are, for example, 95% (or 99.9%, or whatever) sure that the true value lies within this interval.

As a specific example, let us say that we have sampled 100 billion (hopefully) independent and uniformly distributed Monte Carlo realizations of some process, and observed that a certain property is satisfied in 31,895,060,547 of those samples. Each of these samples can be regarded as a Bernoulli trial succeeding with some unknown probability $p$, which we wish to estimate. Thus, the number of successes in $n$ samples is binomially distributed with parameters $n$ and $p$; $n$ is known, and we wish to estimate and characterize the distribution of $p$, given the number of observed successes.

We could now proceed in several ways:

• We could directly apply Bayes' theorem to calculate the conditional distribution of $p$, given the observed number of successes (which will typically be a Beta distribution, if we were smart enough to pick one as our assumed prior distribution of $p$), and then summarize it e.g. by calculating a credible interval for $p$ at whatever confidence level we consider reasonable.

• We could skip the Bayesian calculation, and just directly apply one of the standard formulae for a binomial proportion confidence interval, such as the Wilson or the Jeffreys interval.

• Since we have a lot of samples, and a reasonable proportion of successes, we could simply approximate the distribution of $p$ by a normal distribution having mean equal to the sample mean $\hat p = \frac kn$, where $k$ is the number of observed successes, and variance $\sigma^2 = \frac1n\hat p(1-\hat p)$, and then just quote $\hat p \pm 2\sigma$ as our (slightly conservative) 95% confidence interval for $p$.

Using any of these methods, for the example above, we can obtain a 95% confidence interval of $p \in 0.318951 \pm 2.9\times10^{-6}$ or so. If we wanted to be really sure of getting the true value within the interval, we could widen it to, say $\pm 6.11\sigma = 9.0\times10^{-6}$, and be 99.9999999% sure that the true value of $p$ lies within the interval. (Of course, at those confidence levels, we should start worrying about the non-zero probability of undetected bugs in our Monte Carlo code, or in the random number generator it's using, or in the computer hardware it we ran it on, or in the theory it's based on...)

In fact, for quick and dirty estimation of the accuracy of binomial Monte Carlo results, we can simplify the normal approximation further by noting that $(1-\hat p) \le 1$, and so $n\sigma$ $=$ $\sqrt{n\hat p(1-\hat p)}$ $\le$ $\sqrt{n\hat p}$ $=$ $\sqrt k$. Thus, if we observe a $d$-digit number of successes among our trials, we may roughly estimate that the resulting estimate of the success probability has about $d/2$ significant digits, with the remaining digits being essentially random.

Answer 4

I believe that under certain circumstances, you can do better than just getting a confidence interval. Suppose you decide to estimate $\int_0^1f(x)\,dx$ by $n^{-1}\sum_1^nf(x_i)$ for some random (or pseudorandom) sequence $x_1,x_2,\dots$. If your function $f$ is of bounded variation $V(f)$, then Koksma's Theorem says $$\left|\int f(x)\,dx-{1\over n}\sum_1^nf(x_i)\right|\le V(f)D(x_1,\dots,x_n)$$ where $D$ is the discrepancy of the initial segment of the sequence.