<jupyter_start><jupyter_text>Monte Carlo Estimation: Computing Estimation ErrorWe use the Central Limit Theorem (CLT) to compute how accurate our Monte Carlo estimates are.<jupyter_code>import collections
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as stats
%matplotlib inline<jupyter_output><empty_output><jupyter_text>Example 1We seek the expected value of $\theta$, or the mean of the Gamma distribution defined below.<jupyter_code>m = 10000
a, b = 2, 1/3
theta = stats.gamma.rvs(a, scale=(1/b), size=m)<jupyter_output><empty_output><jupyter_text>We seek the sample mean of theta, which approximately follows a Normal distribution, where the mean of this Normal distribution is the true expected value of $\theta$, and the variance of the Normal distribution is the true variance of $\theta$ divided by m. Thus as we increase m, variance gets smaller. 90% confidence interval for $\theta$Formula for standard error is the standard deviation of the Normal distribution divided by square root of m.The 90% confidence interval for a normal distribution is $\mu$ $\pm$ 2\*standard_error. This is the interval where we are 90% confident that values from the distribution will occur.<jupyter_code>mean = np.mean(theta)
se = np.std(theta) / np.sqrt(m)
print("                   mean: {:.5f}".format(mean))
print("         standard error: {:.5f}".format(se))
ci_90_lb, ci_90_ub = mean - (2 * se), mean + (2 * se)
print("90% confidence interval: ({:.5f}, {:.5f})".format(ci_90_lb, ci_90_ub))<jupyter_output>mean: 5.97243
         standard error: 0.04309
90% confidence interval: (5.88626, 6.05860)<jupyter_text>90% confidence interval for Pr($\theta$ < 5)<jupyter_code>theta_ind = theta < 5
mean = np.mean(theta_ind)
se = np.std(theta_ind) / np.sqrt(m)
print("                   mean: {:.5f}".format(mean))
print("         standard error: {:.5f}".format(se))
ci_90_lb, ci_90_ub = mean - (2 * se), mean + (2 * se)
print("90% confidence interval: ({:.5f}, {:.5f})".format(ci_90_lb, ci_90_ub))<jupyter_output>mean: 0.50480
         standard error: 0.00500
90% confidence interval: (0.49480, 0.51480)<jupyter_text>Error Estimate for Hierarchical ModelWe define our hierarchical model as follow:$$y_i \sim Binomial(10, \phi_i)$$$$\phi_i \sim Beta(2, 2)$$The resulting sample of y can be treated as the marginal distribution. The distribution for (y|$\phi$) is a binomial distribution, but marginal distribution for y is __not binomial__, it is beta-binomial.<jupyter_code>m = int(1e5)
phi = stats.beta.rvs(2, 2, size=m)
y = stats.binom.rvs(10, phi, size=m)
plt.hist(y, density=True)
plt.show()<jupyter_output><empty_output><jupyter_text>Probability of 3 successes, Pr(y=3)<jupyter_code>ind_y = (y == 3)
np.mean(ind_y)<jupyter_output><empty_output><jupyter_text>Marginal Expected Value of y<jupyter_code>np.mean(y)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12281315605318405, 'fraction_numerical': 0.041287613715885234, 'mean_word_length': 5.070063694267516, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '280245', 'license': 'Apache-2.0', 'n_tokens_mistral': 1049, 'n_tokens_neox': 943, 'n_words': 305, 'path': 'techniques-and-models/w01-03-mc-computing-errors.ipynb', 'repo_name': 'sujitpal/bayesian-stats-examples', 'size': '2858'}