017_03.py

sigma2_x = ((X - x_bar)**2).sum()/X.size
sigma2_y = ((Y - y_bar)**2).sum()/Y.size

sigma_x = np.sqrt(sigma2_x)
sigma_y = np.sqrt(sigma2_y)

sigma_xy = (X*Y).sum()/X.size - x_bar*y_bar

r = sigma_xy / (sigma_x*sigma_y)

fig, ax = plt.subplots(facecolor="w")
plt.plot(X, Y, "o", label="observations")
plt.axvline(x_bar, ls="--", label=fr"$\bar{{x}}={x_bar:.2f}$", color="C1")
plt.axhline(y_bar, ls="--", label=fr"$\bar{{y}}={y_bar:.2f}$", color="C2")
X_continuous = np.linspace(min(X), max(X), 100)
sigma_x_1 = (X_continuous>=x_bar-sigma_x)&(X_continuous<=x_bar+sigma_x)
plt.fill_between(
    X_continuous[sigma_x_1],
    min(Y), max(Y),
    color="C1", alpha=.15,
    label=fr"$\sigma_{{x}}={sigma_x:.2f}$"
)
Y_continuous = np.linspace(min(Y), max(Y), 100)
sigma_y_1 = (Y_continuous>=y_bar-sigma_y)&(Y_continuous<=y_bar+sigma_y)
plt.fill_betweenx(
    Y_continuous[sigma_y_1],
    min(X), max(X),
    color="C2", alpha=.15,
    label=fr"$\sigma_{{y}}={sigma_y:.2f}$"
)
plt.text(1.05, .4, f"$r={r:.9f}$", ha="left", va="center", transform=ax.transAxes)
plt.legend(loc="upper left", bbox_to_anchor=(1,1))
plt.show()