Plotting convergence of Monte Carlo estimation of Pi to display convergence rate

Question

The answer to this thread states the variance for 2 Monte Carlo estimators of $\pi$.

Variance of Area and Average Estimators in Monte Carlo Estimation of Pi

The variance for both methods is proportional to $\frac{1}{n}$, so the order of convergence must be $O(\frac{1}{\sqrt{n}})$. I am looking for a visualization of this convergence.

This is what I've done so far; it's a plot of the standard deviation with respect to the sample size. However, I'm not sure if it really captures the $O(\frac{1}{n})$ relationship. Would you have any other suggestion?

# import seaborn as sb
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import math
# sns.set(rc={'text.usetex' : True})


plt.clf()

n = int(1e3)
pi_estimate1 = np.full(n, np.nan)
pi_estimate2 = np.full(n, np.nan)
std1 = np.full(n, np.nan)
std2 = np.full(n, np.nan)
sample_size = 1 + np.arange(n)

for i in range(n):
    np.random.seed(i)
    u1 = np.random.uniform(size = i)
    u2 = np.random.uniform(size = i)
    x = u1**2 + u2**2

    darts = np.zeros(i)
    darts[x <= 1] = 1
    in_circle_proportion = np.mean(darts)
    pi_estimate1[i] = 4 * in_circle_proportion
    std1[i] = 4*np.std(darts)/np.sqrt(i)

    sample = 4 * np.sqrt(1 - u1 ** 2)
    pi_estimate2[i] = np.mean(sample)
    std2[i] = np.std(sample)/np.sqrt(i)


plt.scatter(sample_size, std1, s = 3, label = 'Area method')
plt.scatter(sample_size, std2, s = 3, label = 'Monte Carlo Integration')
plt.xlabel("Sample Size (n)")
plt.ylabel('Standard Deviation')
plt.ylim(top = 1)
plt.ylim(bottom = 0)
plt.xlim(left = 0)
plt.xlim(right = 200)
plt.legend()
plt.show()

log of the standard deviation or log of the sample size? and why should I do this? — Iterator516, Jan 22 '20 at 20:00
You are expecting a relationship of the form $d_n \le C {1 \over \sqrt{n}}$, so if you plot $\log d_k $ vs. $\log n$ you would expect an upper bound that is a line of slope $-1$ (since $\log d_k \le \log C - {1 \over 2} (\log n)$). — copper.hat, Jan 22 '20 at 20:08
I plotted log(std) vs. log(n) as you suggested, and I see that the slope is roughly -0.5, which is what you showed: the slope of -1/2 for log(n). Is this what you mean? — Iterator516, Jan 22 '20 at 20:15
Yes. That is a typical way of demonstrating rate of convergence. — copper.hat, Jan 22 '20 at 20:15
Thank you. I provided a solution with the code and the plot. Does this look right to you? — Iterator516, Jan 22 '20 at 20:30
Looks good to me. Except I would use $\log$ axes or label the axes with $\log n$ and simiiar. — copper.hat, Jan 22 '20 at 20:39

Iterator516 · Answer 1 · 2020-01-22T21:08:21.893

As @copper.hat suggested, I plotted the log of the standard deviation with the log of the sample size.

# import seaborn as sb
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import math
# sns.set(rc={'text.usetex' : True})


plt.clf()

n = int(1e3)
pi_estimate1 = np.full(n, np.nan)
pi_estimate2 = np.full(n, np.nan)
std1 = np.full(n, np.nan)
std2 = np.full(n, np.nan)
sample_size = 1 + np.arange(n)

for i in range(n):
    np.random.seed(i)
    u1 = np.random.uniform(size = i)
    u2 = np.random.uniform(size = i)
    x = u1**2 + u2**2

    darts = np.zeros(i)
    darts[x <= 1] = 1
    in_circle_proportion = np.mean(darts)
    pi_estimate1[i] = 4 * in_circle_proportion
    std1[i] = 4*np.std(darts)/np.sqrt(i)

sample = 4 * np.sqrt(1 - u1 ** 2)
pi_estimate2[i] = np.mean(sample)
std2[i] = np.std(sample)/np.sqrt(i)


plt.scatter(np.log(sample_size), np.log(std1), s = 30, label = 'Area method')
plt.scatter(np.log(sample_size), np.log(std2), s = 30, label = 'Monte Carlo Integration')
plt.xlabel("ln(Sample Size)")
plt.ylabel('ln(Standard Deviation)')
plt.legend(fontsize = 'large')
plt.show()

Plotting convergence of Monte Carlo estimation of Pi to display convergence rate

1 Answers1