1

In a computational physics course, I was asked to do direct-sampling for the numerical value of $\pi$

and then I estimated the standard deviation of $\pi$ , math.sqrt(var / n_trials) further by running Python. The standard deviation however is $\longrightarrow 1.64$ as the total $N$ (n_trial) increases.

I do not quite understand the meaning of this $1.64$, or more generally the standard deviation in this kind of algorithm. Does that mean the smaller the standard deviation is, the "better" this algorithm is? Or it has no special meaning, as long as it converges, I am having a reasonable method?

Here is the code for Python provided by the course:

import random, math
n_trials = 400000                                #n_trial is N
n_hits = 0
var = 0.0
for iter in range(n_trials):
    x, y = random.uniform(-1.0, 1.0), random.uniform(-1.0, 1.0)
    Obs = 0.0
    if x**2 + y**2 < 1.0:
        n_hits += 1
        Obs = 4.0
    var += (Obs - math.pi)**2
print(4.0 * n_hits / float(n_trials), math.sqrt(var / n_trials)) #estimated pi and standard deviation respectively

p.s. the math.pi is just for convenience, we will have the standard deviation = 1.64 too if we use mean value of the observables instead.

Shing
  • 113
  • 5

1 Answers1

1

Obs (call it $X$) is a RV that takes values $4$ and $0$ with probabilities $\pi/4$ and $1-\pi/4$ respectively. Expected value of it is $\pi$. Variance of it is $$E[X^2]=16\times\pi/4=4\pi\rightarrow \operatorname{var}(X)=4\pi-\pi^2=\pi(4-\pi)$$ And, the deviation is $\sigma(X)=\sqrt{\pi(4-\pi)}\approx 1.64$. I'm not sure why you calculate this deviation but it is certainly not related to your algorithm being good or bad.

A might be better judge of the algorithm will be the variance/deviation of your $\pi$ estimate, i.e. $\hat \pi=4\times \frac{\text{hits}}{\text{trials}}$. Here hits is a binomial RV (call it $Y$) with parameters n=n_trials and $p=\pi/4$. Binomial RV's have variance $np(1-p)$, therefore we have: $$\operatorname{var}(\hat\pi)=\operatorname{var}\left(4Y/n\right)=\frac{16}{n^2}(np(1-p))=\frac{4}{n}\pi(1-\pi/4)\approx (1.64)^2/n$$ which indeed decreases with increasing $n$.

gunes
  • 49,700
  • 3
  • 39
  • 75
  • Thanks for the answer. sorry, my bad. I did not mean I calculated the SD, but actually I got the SD ~ 1.64 via running the Python codes || Is `n=trials` --> `n_trials` ? – Shing Jan 18 '20 at 19:08
  • Yes, I changed it. – gunes Jan 18 '20 at 19:24
  • I have upvoted, but I am still digesting your answer, once I have learned your answer completely, I will accept your answer for sure! btw, would you mind if I edit your notation such that your math/code will be in consistence with my question? (I can edit my question instead, if you prefer) – Shing Jan 19 '20 at 11:32
  • So if I have a "steeper" $\text{var}(\hat \pi)$, then I have a better algorithm? – Shing Feb 07 '20 at 17:05
  • It's a good sign. In general, you can't only look at the variance, e.g. $\hat \pi = 3$ has variance zero but a poor estimate. However, if the estimator is also unbiased, it's a better sign. Variance going zero provides consistency. So, steepness can be thought as a measure of consistency as well. You can read more in https://stats.stackexchange.com/questions/418417/why-do-we-need-an-estimator-to-be-consistent/418422#418422 – gunes Feb 07 '20 at 17:19