Monte Carlo approximation of pi using normally distributed points

Question

I'm trying to calculate pi using the classic Monte Carlo method. But, instead of using uniformly distributed points, I want to use a normal distribution with mean value centred on the origin ($\mu = 0$) and standard deviation of $\sigma = 0.6$.

If the darts are thrown uniformly in the square from $(-1,-1)$ to $(1,1)$, this distance is less than $1$ with a probability of $P(\textrm{inside circle}) = \frac{\pi}{4}$. How would this change for a normal distribution? How would you calculate $\pi$?

My attempt:

For the normal distribution described above $P(\textrm{inside circle}) = 0.905$ (using the look-up table for $Z$). But this expression doesn't include $\pi$, so I'm not sure where to go from here.

Any help would be much appreciated.

Answer 1

You can use distribution function transformation (cf Casella Berger or any standard statistics texts) to get uniform distribution from normal. Basically, with continuous r.v. X, we have $F_X(X) \sim uniform (0,1)$, where $ F_X $ is distribution function of $X$. Here, you need to first estimate the distribution function using the empirical distribution function $\hat F_X(x)=\sum_{i=1}^n1\{X_i\le x\}/n$. Glivenko-Cantelli theorem guarantee good properties of this estimate.

Answer 2

tl;dr

yet another way to approximate pi with normally distributed points. See example at the end.

Preamble

I recently stumbled across a very similar problem. I am not sure your proposed approach (sampling the unit square) is feasible, because your points are not uniformly distributed and you get more mass in the inscribed circle. Hereby is a valid alternative, which is still based on a Monte Carlo approximation. (Also - needless to say - Chuck's approach is totally valid, the below is just different method which I found quite interesting and - if anything - it saves you from the calculation of the empirical CDF)

I think we can all agree that using a $Uniform(0,1)$ distribution is more intuitive, however, if you really wanted to unnecessarily flex and stick to a $N(\mu,\sigma)$, you can exploit the fact that the pdf of a $N(\mu,\sigma)$ is a function of $\pi$. So the problem now becomes a matter of inverting the pdf equation to get yourself what you want.

Methodology

Below is the pdf of a $N(\mu,\sigma)$

$$p(x) = \frac{1}{\sigma \sqrt{2\pi}}e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2} \tag{1}$$

You can think of $(1)$ as of 2 parts: $\frac{1}{\sqrt{2\pi}}$ and $\frac{1}{\sigma}e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2}$

Although the pdf of a $N(\mu,\sigma)$ does not integrate analytically (because $e^{-x^2} $ does not), its first central moment does (you can try either by parts or by substitution. It is easier to get an intuition of that if you start with the simpler case $\mu=0$ and $\sigma=1$)

$$ \frac{1}{\sqrt{2\pi}} \int (x-\mu) \cdot \frac{1}{\sigma}e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2} \, dx = \frac{1}{\sqrt{2\pi}} \cdot \sigma\,\big(-e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2}\,\big)+K \tag{2} $$

You can then write the following equation:

$$ \frac{1}{\sqrt{2\pi}} \cdot \frac{A}{B} = \frac{C}{B} \tag{3} $$

Where:

• $A := \int_a^b x \cdot \frac{1}{\sigma}e^{-\frac{1}{2}(\frac{x-\mu}{\sigma})^2} \, dx $
• $B := \int_a^b p(x) \, dx$ i.e. the probability mass of between $a$ and $b$
• $C := \int_a^b x \cdot p(x) \, dx$ i.e. the first raw moment of the pdf between $a$ and $b$

Once you realize the above and you observe that you can numerically compute

• $\frac{C}{B}$ i.e. the normalized central moment (it is the sample average less $\mu$ after all)
• $B$ i.e. the probability mass (it is the sample frequency)

of your $N(\mu,\sigma)$ between any two points $a$ and $b$.

You then simply rearrange $(3)$ to have

$$ \pi = \frac{\big(\frac{A}{B \, C/B}\big)^2}{2} \tag{4}$$

Keep in mind that you will have to choose your extremes of integration wisely (most notably, you will definitely have to make sure $a\not=-b$ otherwise $(4)$ is undefined). Quantity $B$ is required to normalize A and C (remember that you are integrating between $a$ and $b$)

If the value B in $(4)$ seems redundant to you, remember that you cannot compute C directly (i.e. with its analytical form as per $(2)$, because you do not know the value of $\pi$)

Implementation

At this GitHub link, my Python implementation of the above. Try it (e.g. with your parameters). You will quickly realize that this is not the most efficient way to accurately compute $\pi$ to any sort of precision

Example

With the values from the original question: $\mu = 0$, $\sigma = 0.6$ (and assuming you choose $a=0$, $b=\infty$) you get

• $A = \int_0^\infty x \cdot \frac{1}{0.6}e^{-\frac{1}{2}(\frac{x}{0.6})^2} \, dx = 0.6 $ by the analitical solution in $(2)$
• $B = \int_0^\infty \frac{1}{0.6\sqrt{2\pi}}e^{-\frac{1}{2}(\frac{x}{0.6})^2} \, dx = 0.5$ trivial; otherwise can be calculated numerically (i.e. numbers of sampled points between $a$ and $b$ over the total number of sampled )
• $C/B = {\displaystyle \int_0^\infty} \frac{ x \cdot \frac{1}{0.6\sqrt{2\pi}}e^{-\frac{1}{2}(\frac{x}{0.6})^2}}{\int_0^\infty \frac{1}{0.6\sqrt{2\pi}}e^{-\frac{1}{2}(\frac{x}{0.6})^2}\, dx} \, dx \approx 0.47873$ calculated numerically (i.e. sample average between $a$ and $b$)

Replacing in $(4)$ you have that:

$$\pi \approx \frac{\big(\frac{0.6}{0.5 \cdot 0.47873 }\big)^2}{2} = 3.14160$$

If you really wanted to use normally scattered point (as per your example) in the XY plane as your PseudoRandomNumberGenerator engine, imagine to take the X (or Y) component of your coordinates only.