<jupyter_start><jupyter_text>Code boxes Chapter 3<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
import pymc3 as pm<jupyter_output><empty_output><jupyter_text>Py code 3.1<jupyter_code>Pr_Positive_Vampire = 0.95
Pr_Positive_Mortal = 0.01
Pr_Vampire = 0.001
Pr_Positive = Pr_Positive_Vampire * Pr_Vampire + Pr_Positive_Mortal * (1 - Pr_Vampire)

Pr_Vampire_Positive = Pr_Positive_Vampire * Pr_Vampire / Pr_Positive
print(Pr_Vampire_Positive)<jupyter_output>0.08683729433272395<jupyter_text>Py code 3.2<jupyter_code>p_grid = np.linspace(0, 1, num=1000)
prob_p = np.ones(1000)
prob_data = stats.binom.pmf(k=6, n=9, p=p_grid)
posterior = prob_data * prob_p
posterior /= np.sum(posterior)<jupyter_output><empty_output><jupyter_text>Py code 3.3<jupyter_code>samples = np.random.choice(p_grid, p=posterior, size=10000, replace=True)<jupyter_output><empty_output><jupyter_text>Py code 3.4The python code below actually reproduces a figure similar to the one shown in the book, Fig 3.1. The code however would make a more simple figure.<jupyter_code>plt.plot(samples, 'o', alpha=0.2)
plt.xlabel('sample number')
plt.ylabel('proportion water (p)')
plt.ylim([0, 1])<jupyter_output><empty_output><jupyter_text>Py code 3.5Instead of the `rethinking` package we are again using the kernel density plot from `seaborn`.<jupyter_code>sns.kdeplot(samples)
plt.xlabel('proportion water (p)')
plt.ylabel('Density')<jupyter_output><empty_output><jupyter_text>Py code 3.6<jupyter_code>np.sum(posterior[np.where(p_grid < 0.5)])<jupyter_output><empty_output><jupyter_text>Py code 3.7<jupyter_code>np.sum(samples<0.5) / 1e4<jupyter_output><empty_output><jupyter_text>*Note*: In Py code 3.6 we summed up the values of posterior at the indices where p_grid has values < 0.5. Here we count the number of samples that are < 0.5. `samples<0.5` on the numpy array `samples` simply creates a boolean mask, building the sum over it gives us the number of samples because `True` is interpreted as 1 and `False` is interpreted as `0`. This is the same as in `R`, see the "overthinking" paragraph below these examples in the book.  Py code 3.8<jupyter_code>np.sum(np.logical_and(samples > 0.5, samples < 0.75)) / 1e4<jupyter_output><empty_output><jupyter_text>Py code 3.9<jupyter_code>np.quantile(samples, 0.8)<jupyter_output><empty_output><jupyter_text>Py code 3.10<jupyter_code>np.quantile(samples, [0.1, 0.9])<jupyter_output><empty_output><jupyter_text>Py code 3.11<jupyter_code>p_grid = np.linspace(0, 1, num=1000)
prior = np.ones(1000)
likelihood = stats.binom.pmf(k=3, n=3, p=p_grid)
posterior = likelihood * prior
posterior /= np.sum(posterior)

samples = np.random.choice(p_grid, p=posterior, size=10000, replace=True)<jupyter_output><empty_output><jupyter_text>Py code 3.12<jupyter_code>np.quantile(samples, [0.25, 0.75])<jupyter_output><empty_output><jupyter_text>*Note*: Using the quantile function is identical to using the percentile function as following:<jupyter_code>np.percentile(samples, [25, 75])<jupyter_output><empty_output><jupyter_text>Py code 3.13<jupyter_code>pm.stats.hpd(samples, credible_interval=0.5)<jupyter_output><empty_output><jupyter_text>Py code 3.14<jupyter_code>p_grid[np.argmax(posterior)]<jupyter_output><empty_output><jupyter_text>Py code 3.15<jupyter_code>stats.mode(samples)[0]<jupyter_output><empty_output><jupyter_text>Py code 3.16<jupyter_code>np.mean(samples)
np.median(samples)<jupyter_output><empty_output><jupyter_text>Py code 3.17<jupyter_code>np.sum(posterior * np.abs(0.5 - p_grid))<jupyter_output><empty_output><jupyter_text>Py code 3.18<jupyter_code>loss = [np.sum(posterior * np.abs(x - p_grid)) for x in p_grid]<jupyter_output><empty_output><jupyter_text>Py code 3.19<jupyter_code>p_grid[np.argmin(loss)]<jupyter_output><empty_output><jupyter_text>Py code 3.20<jupyter_code>stats.binom.pmf(k=range(3), n=2, p=0.7)<jupyter_output><empty_output><jupyter_text>Py code 3.21<jupyter_code>stats.binom.rvs(n=2, p=0.7, size=1)<jupyter_output><empty_output><jupyter_text>Py code 3.22<jupyter_code>stats.binom.rvs(n=2, p=0.7, size=10)<jupyter_output><empty_output><jupyter_text>Py code 3.23<jupyter_code>dummy_w = stats.binom.rvs(n=2, p=0.7, size=int(1e5))
[(dummy_w == i).mean() for i in range(3)]<jupyter_output><empty_output><jupyter_text>*Note*: `dummy_w == i` makes a mask on the numpy array. as before, all `True` values are interpreted as ones, all `False` as zeros when calculating the mean.  Py code 3.24<jupyter_code>dummy_w = stats.binom.rvs(n=9, p=0.7, size=int(1e5))
plt.hist(dummy_w, bins=np.arange(-0.5, 10, 1), rwidth=0.1)
plt.xlabel('dummy water count')
plt.ylabel('Frequency')<jupyter_output><empty_output><jupyter_text>*Note*: Here we are forcing 10 bins from 0 to 9 that are right on the integer by giving a 11 bin ranges. See `plt.hist?` for help. Alternatively we could also specify 50 - 100 bins, that way we would likely have a bin that is very close to an integer value.  Py code 3.25<jupyter_code>w = stats.binom.rvs(n=9, p=0.6, size=int(1e4))

# Also plot the distribution, as above:
plt.hist(w, bins=50)
plt.show()<jupyter_output><empty_output><jupyter_text>Py code 3.26 PreparationThe first cell below reruns the samples from before. The current `samples` are the ones from Figure 3.3 and are not the ones with `n=9`, `k=6`.<jupyter_code>p_grid = np.linspace(0, 1, num=1000)
prob_p = np.ones(1000)
prob_data = stats.binom.pmf(k=6, n=9, p=p_grid)
posterior = prob_data * prob_p
posterior /= np.sum(posterior)
samples = np.random.choice(p_grid, p=posterior, size=10000, replace=True)<jupyter_output><empty_output><jupyter_text>Code<jupyter_code>w = stats.binom.rvs(n=9, p=samples, size=int(1e4))

# also plot
plt.hist(w, bins=50)
plt.show()<jupyter_output><empty_output><jupyter_text>Py code 3.27<jupyter_code>p_grid = np.linspace(0, 1, num=1000)
prior = np.ones(1000)
likelihood = stats.binom.pmf(k=6, n=9, p=p_grid)
posterior = likelihood * prior
posterior /= np.sum(posterior)
np.random.seed(100)
samples = np.random.choice(p_grid, p=posterior, size=int(1e4), replace=True)<jupyter_output><empty_output><jupyter_text>Py code 3.28<jupyter_code>birth1 = np.array([1,0,0,0,1,1,0,1,0,1,0,0,1,1,0,1,1,0,0,0,1,0,0,0,1,0,
                   0,0,0,1,1,1,0,1,0,1,1,1,0,1,0,1,1,0,1,0,0,1,1,0,1,0,
                   0,0,0,0,0,0,1,1,0,1,0,0,1,0,0,0,1,0,0,1,1,1,1,0,1,0,
                   1,1,1,1,1,0,0,1,0,1,1,0,1,0,1,1,1,0,1,1,1,1])
birth2 = np.array([0,1,0,1,0,1,1,1,0,0,1,1,1,1,1,0,0,1,1,1,0,0,1,1,1,0,
                   1,1,1,0,1,1,1,0,1,0,0,1,1,1,1,0,0,1,0,1,1,1,1,1,1,1,
                   1,1,1,1,1,1,1,1,1,0,1,1,0,1,1,0,1,1,1,0,0,0,0,0,0,1,
                   0,0,0,1,1,0,0,1,0,0,1,1,0,0,0,1,1,1,0,0,0,0])<jupyter_output><empty_output><jupyter_text>Py code 3.29Not applicable, the two numpy arrays can be copied from above.  Py code 3.30<jupyter_code>np.sum(birth1) + np.sum(birth2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15275567770866483, 'fraction_numerical': 0.07117025893244612, 'mean_word_length': 6.892694063926941, 'pattern_counts': {'":': 0, '<': 134, '<?xml version=': 0, '>': 129, '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': '32', 'dir': 'jupyter-structured-clean-dedup', 'id': '653719', 'license': 'MIT', 'n_tokens_mistral': 3206, 'n_tokens_neox': 2959, 'n_words': 518, 'path': 'code_boxes/Chapter_3.ipynb', 'repo_name': 'trappitsch/PyThinking', 'size': '6913'}