#389 Calculating pi with Monte Carlo Simulation

# Calculating pi with Monte Carlo Simulation
import random
import numpy as np
import matplotlib.pyplot as plt

random.seed(2021)
pi_values = list()
num_experiments = 1000
num_rounds = 20
num_points = 1000
edge = 10
for r in range(num_rounds):
    for p in range(num_experiments):
        in_circle = 0
        for g in range(num_points):
            x, y = (random.random() - 0.5)*edge, (random.random() - 0.5)*edge
            if x**2 + y**2 <= (edge/2)**2:
                in_circle = in_circle + 1
        pi = in_circle/num_points * 4
        pi_values.append(pi)
print(np.mean(pi_values),np.pi)        

random.seed(2021)
edge = 10
num_points = 1000
with plt.xkcd():
    fig, ax = plt.subplots(figsize=(6, 6))
# plt.axis("off")
    plt.axis("equal")
    ax.set_xlim(-edge/2, edge/2)
    ax.set_ylim(-edge/2, edge/2)
    xs_in, ys_in = list(), list()
    xs_out, ys_out = list(), list()
    for g in range(num_points):
        x, y = (random.random() - 0.5)*edge, (random.random() - 0.5)*edge
        if x**2 + y**2 <= (edge/2)**2:
            xs_in.append(x)
            ys_in.append(y)
        else:
            xs_out.append(x)
            ys_out.append(y)
    ax.scatter(xs_in, ys_in, color="r")
    ax.scatter(xs_out, ys_out, color="b")
    circle = plt.Circle((0, 0), edge/2, fill=False, color="g", lw=3)
    ax.add_patch(circle)
    ax.set_title("An experiment with 1000 points", fontsize=20)
    plt.savefig("ex389_1.png", dpi=100)
    plt.show()
    
    
with plt.xkcd():
    fig, ax = plt.subplots(figsize=(12, 6))
    ax.hist(pi_values, bins=50, rwidth=0.8)
    pi_mean = np.mean(pi_values)
    pi_median = np.median(pi_values)
    # pi_mean and pi_median are very close.
    
    pi_std = np.std(pi_values)
    pi_quartiles = np.quantile(pi_values, [0, 0.25, 0.5, 0.75, 1])
    ax.set_title(
        "Statistics of the Experiments", fontsize=20)
    line_1 = ax.axvline(pi_mean, color='red', lw=1)
    line_2 = ax.axvline(pi_quartiles[1],
                        color='purple',
                        lw=3,
                        linestyle="dotted")
    line_3 = ax.axvline(pi_quartiles[3],
                        color='green',
                        lw=3,
                        linestyle="dashed")
    std_bar = ax.errorbar(pi_mean, 1200,
                          xerr=pi_std,
                          capsize=5,
                          elinewidth=3,
                          markeredgewidth=2,
                          linestyle=":")
    ax.legend([line_1, line_2, line_3, std_bar],
              ["median = {}".format(pi_median),
               "first quartile = {}".format(pi_quartiles[1]),
               "third quartile = {}".format(pi_quartiles[3]),
               "standard deviation =\n {}".format(round(pi_std, 4))],
              fontsize=18)    
    plt.savefig("ex389_2.png", dpi=100)
    plt.show()