<jupyter_start><jupyter_text>Chapter 3. Sampling the Imaginary<jupyter_code>import os

import arviz as az
import matplotlib.pyplot as plt
import pandas as pd
from scipy.stats import gaussian_kde

import jax.numpy as jnp
from jax import random, vmap

import numpyro
import numpyro.distributions as dist

if "SVG" in os.environ:
    %config InlineBackend.figure_formats = ["svg"]
az.style.use("arviz-darkgrid")<jupyter_output><empty_output><jupyter_text>Code 3.1<jupyter_code>Pr_Positive_Vampire = 0.95
Pr_Positive_Mortal = 0.01
Pr_Vampire = 0.001
tmp = Pr_Positive_Vampire * Pr_Vampire
Pr_Positive = tmp + Pr_Positive_Mortal * (1 - Pr_Vampire)
Pr_Vampire_Positive = tmp / Pr_Positive
Pr_Vampire_Positive<jupyter_output><empty_output><jupyter_text>Code 3.2<jupyter_code>p_grid = jnp.linspace(start=0, stop=1, num=1000)
prob_p = jnp.repeat(1, 1000)
# Throw the world 9 times, 6 times lands on water, what is the amount of water on the world?
prob_data = jnp.exp(dist.Binomial(total_count=9, probs=p_grid).log_prob(6))
posterior = prob_data * prob_p
posterior = posterior / jnp.sum(posterior)<jupyter_output><empty_output><jupyter_text>Code 3.3<jupyter_code>plt.scatter(p_grid, posterior)
# Categorical samples from 0,...,1 with probabilities of the posterior
samples = p_grid[dist.Categorical(probs=posterior).sample(random.PRNGKey(0), (10000,))]<jupyter_output><empty_output><jupyter_text>Code 3.4<jupyter_code>plt.scatter(range(len(samples)), samples, alpha=0.2)
plt.show()<jupyter_output><empty_output><jupyter_text>Code 3.5<jupyter_code>az.plot_density({"": samples}, hdi_prob=1)
plt.show()<jupyter_output><empty_output><jupyter_text>Code 3.6<jupyter_code># add up posterior probability where p < 0.5
jnp.sum(posterior[p_grid < 0.5])<jupyter_output><empty_output><jupyter_text>Code 3.7<jupyter_code>jnp.sum(samples < 0.5) / 1e4<jupyter_output><empty_output><jupyter_text>Code 3.8<jupyter_code>jnp.sum((samples > 0.5) & (samples < 0.75)) / 1e4<jupyter_output><empty_output><jupyter_text>Code 3.9<jupyter_code>jnp.quantile(samples, 0.8)<jupyter_output><empty_output><jupyter_text>Code 3.10<jupyter_code>jnp.quantile(samples, jnp.array([0.1, 0.9]))<jupyter_output><empty_output><jupyter_text>Code 3.11<jupyter_code>p_grid = jnp.linspace(start=0, stop=1, num=1000)
prior = jnp.repeat(1, 1000)
likelihood = jnp.exp(dist.Binomial(total_count=3, probs=p_grid).log_prob(3))
posterior = likelihood * prior
posterior = posterior / jnp.sum(posterior)
samples2 = p_grid[dist.Categorical(probs=posterior).sample(random.PRNGKey(0), (10000,))]<jupyter_output><empty_output><jupyter_text>Code 3.12<jupyter_code>jnp.percentile(samples2, q=(25, 75))<jupyter_output><empty_output><jupyter_text>Code 3.13<jupyter_code>numpyro.diagnostics.hpdi(samples2, prob=0.5)<jupyter_output><empty_output><jupyter_text>Code 3.14<jupyter_code>p_grid[jnp.argmax(posterior)]<jupyter_output><empty_output><jupyter_text>Code 3.15<jupyter_code>samples2[jnp.argmax(gaussian_kde(samples2, bw_method=0.01)(samples))]<jupyter_output><empty_output><jupyter_text>Code 3.16<jupyter_code>print(jnp.mean(samples2))
print(jnp.median(samples2))<jupyter_output>0.8011086
0.8428428<jupyter_text>Code 3.17<jupyter_code>jnp.sum(posterior * jnp.abs(0.5 - p_grid))<jupyter_output><empty_output><jupyter_text>Code 3.18<jupyter_code>loss = vmap(lambda d: jnp.sum(posterior * jnp.abs(d - p_grid)))(p_grid)
plt.scatter(p_grid, loss);<jupyter_output><empty_output><jupyter_text>Code 3.19<jupyter_code>p_grid[jnp.argmin(loss)]<jupyter_output><empty_output><jupyter_text>Code 3.20<jupyter_code>jnp.exp(dist.Binomial(total_count=2, probs=0.7).log_prob(jnp.arange(3)))<jupyter_output><empty_output><jupyter_text>Code 3.21<jupyter_code>dist.Binomial(total_count=2, probs=0.7).sample(random.PRNGKey(0))<jupyter_output><empty_output><jupyter_text>Code 3.22<jupyter_code>dist.Binomial(total_count=2, probs=0.7).sample(random.PRNGKey(2), (10,))<jupyter_output><empty_output><jupyter_text>Code 3.23<jupyter_code>dummy_w = dist.Binomial(total_count=2, probs=0.7).sample(random.PRNGKey(0), (100000,))
jnp.unique(dummy_w, return_counts=True)[1] / 1e5<jupyter_output><empty_output><jupyter_text>Code 3.24<jupyter_code>dummy_w = dist.Binomial(total_count=9, probs=0.7).sample(random.PRNGKey(0), (100000,))
ax = az.plot_dist(dummy_w.copy(), kind="hist", hist_kwargs={"rwidth": 0.1})
ax.set_xlabel("dummy water count", fontsize=14)
plt.show()<jupyter_output><empty_output><jupyter_text>Code 3.25<jupyter_code>w = dist.Binomial(total_count=9, probs=0.6).sample(random.PRNGKey(0), (int(1e4),))
plt.hist(w)<jupyter_output><empty_output><jupyter_text>Code 3.26<jupyter_code>plt.hist(samples)
w = dist.Binomial(total_count=9, probs=samples).sample(random.PRNGKey(0))
plt.hist(w)<jupyter_output><empty_output><jupyter_text>Practice   Code 3.27<jupyter_code>p_grid = jnp.linspace(start=0, stop=1, num=1000)
prior = jnp.repeat(1, 1000)
likelihood = jnp.exp(dist.Binomial(total_count=9, probs=p_grid).log_prob(6))
posterior = likelihood * prior
posterior = posterior / jnp.sum(posterior)
samples = p_grid[dist.Categorical(posterior).sample(random.PRNGKey(100), (10000,))]
#3E1
sum(samples < 0.2) / len(samples)<jupyter_output><empty_output><jupyter_text>Code 3.28<jupyter_code># fmt: off
birth1 = [
    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 = [
    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>Code 3.29<jupyter_code>homeworkch3 = pd.read_csv("../data/homeworkch3.csv")<jupyter_output><empty_output><jupyter_text>Code 3.30<jupyter_code>sum(birth1) + sum(birth2)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16836065573770492, 'fraction_numerical': 0.07754098360655738, 'mean_word_length': 7.346101231190151, 'pattern_counts': {'":': 2, '<': 129, '<?xml version=': 0, '>': 125, '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': '31', 'dir': 'jupyter-structured-clean-dedup', 'id': '215297', 'license': 'MIT', 'n_tokens_mistral': 3128, 'n_tokens_neox': 2743, 'n_words': 457, 'path': 'notebooks/03_sampling_the_imaginary.ipynb', 'repo_name': 'pim-hoeven/rethinking-numpyro', 'size': '6100'}