<jupyter_start><jupyter_text>COMP 135 day03 Lab: Beta distributions and the Beta-Bernoulli model for coin tosses Outline* **Part 1: Beta Distribution visualization*** * Examine plot of Beta PDF at various parameters a>0, b>0.* * Exercises 1a-1c: Inspect samples from several possible Beta distributions. Make connections to the PDF plots.* **Part 2: Beta-Bernoulli Posteriors*** * Exercise 2a: Write function to determine the posterior parameters after seeing an array of binary obervations* * Exercise 2b: Show the posterior after seeing 1, 10, 100, 1000, 10000 examples of a coin* * Exercise 2d: Show the posterior after seeing 1, 10, 100, 1000, 10000 examples, with a different (misinformed) prior Takeaways* Scipy.stats has *many* great resources for evaluating PDFs and for sampling* * Evaluate pdf with `scipy.stats.(hyperparameters).pdf(values)` (can use the `logpdf` instead of `pdf`)* * Draw samples with `scipy.stats.(hyperparameters).rvs(size=(nrows,ncols), random_state=SEED)`* Visualization is an excellent way to understand the behavior of distributions* Analyzing "fake data" drawn from the model is a great debugging strategy.<jupyter_code>import numpy as np
import scipy.stats
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style("whitegrid")
sns.set_context("notebook", font_scale=1.25)<jupyter_output><empty_output><jupyter_text>Part 1: Beta Distribution Visualized Beta distribution as a PDF  Formal definition:  $$\text{BetaPDF}(\mu | a, b) = c(a, b) \cdot \mu^{a-1} (1-\mu)^{b-1}$$  Standard Numerical implementationsThankfully, you can use the built-in functions from SciPy to compute this function:* `stats.beta.pdf(\cdot)` computes the pdf$* `stats.beta.logpdf(\cdot)` computes the natural log of the pdf* * Will often be more numerically stable than doing `np.log(stats.beta.pdf(\cdot))`See   Visualization: PDF across different values of hyperparameters a, bWe'll make a grid plot, where each grid cell shows the PDF for a single $a, b$ pair.<jupyter_code>agrid_G = np.asarray([0.01, 0.1, 0.5, 1.0, 2.0, 10.0, 99.0])
bgrid_H = np.asarray([0.01, 0.1, 0.5, 1.0, 2.0, 10.0, 99.0])
G = agrid_G.size
H = bgrid_H.size

mu_vals = np.linspace(0.0001, 0.9999, 101)

fig, axgrid = plt.subplots(nrows=len(agrid_G), ncols=len(bgrid_H), sharex=True, sharey=True, figsize=(12, 12));
axgrid = axgrid[::-1] # flip order of rows

for gg in range(G):
    for hh in range(H):
        # Evaluate Beta PDF at each candidate value of mu in "mu_vals"
        axgrid[gg, hh].plot(mu_vals, scipy.stats.beta(agrid_G[gg], bgrid_H[hh]).pdf(mu_vals), '-')
        if hh == 0:
            axgrid[gg, hh].set_ylabel('a = %.2g\n\n' % agrid_G[gg]);
        if gg == G - 1:
            axgrid[gg, hh].set_title('b = %.2g' % bgrid_H[hh]);
        if gg == 0:
            axgrid[gg, hh].set_xlabel('$\mu$');

        if agrid_G[gg] == 1.0:
            axgrid[gg, hh].spines['left'].set_linewidth(3)
            axgrid[gg, hh].spines['left'].set_color('darkgray')
        if bgrid_H[hh] == 1.0:
            axgrid[gg, hh].spines['bottom'].set_linewidth(3)
            axgrid[gg, hh].spines['bottom'].set_color('darkgray')

H = 3.0
B = 0.05
sns.despine(fig)
axgrid[0,0].set_ylim(0, H); axgrid[0,0].set_yticks(np.arange(0, H, 1.0));
axgrid[0,0].set_xlim(-B, 1+ B); axgrid[0,0].set_xticks([0, 0.5, 1.0]); axgrid[0,0].set_xticklabels(['0', '0.5', '1']);

#plt.savefig('day03-beta_pdf_viz.pdf', bbox_inches='tight', pad_inches=0);<jupyter_output><empty_output><jupyter_text>Discussion 1a: Inspect 20 samples from a Beta(10, 10). What do you notice about them?<jupyter_code>a = 10
b = 10
mu_samples = scipy.stats.beta(a,b).rvs(size=20, random_state=0)

np.set_printoptions(precision=3, suppress=True)
print(mu_samples)<jupyter_output>[0.597 0.707 0.459 0.549 0.509 0.625 0.238 0.762 0.504 0.603 0.502 0.534
 0.494 0.377 0.395 0.464 0.477 0.584 0.327 0.546]<jupyter_text>TODO what do you notice? How does what you see agree or disagree with the PDF visualizations we created above?  Exercise 1b: Draw 20 samples from a Beta(0.1, 0.1). What do you notice? TODO write code to draw samples, and then discuss  Exercise 1c: Draw 20 samples from a Beta (0.001, 0.001). What do you notice? TODO write code to draw samples, and then discuss  Part 2: Posterior visualization in the Beta-Bernoulli model We'll assume the model discussed in class\begin{align}p(\mu) &= \text{Beta}( \mu | a, b)\\p(x_1, \ldots x_N | \mu) &= \prod_{n=1}^N \text{BernPMF}( x_n | \mu)\end{align}  Exercise 2a: Calculating posterior parametersUsing the formulas discussed in the textbook / video lecture, complete the function below to determine the parameters of the Beta posterior distribution after seeing N coin flip binary outcomes in array `x_N`<jupyter_code>def calc_posterior_parameters(x_N, a, b):
    ''' Compute posterior parameters
    
    Args
    ----
    x_N : 1D array, size N
        Contains binary observations for each of N coin tosses
    a : float, positive
        Prior parameter, number of heads "pseudocounts"
    b : float, positive
        Prior parameter, number of tails "pseudocounts"
    
    Returns
    -------
    ahat : positive scalar
    bhat : positive scalar
    '''
    ahat = a # TODO fixme
    bhat = b # TODO fixme
    return ahat, bhat<jupyter_output><empty_output><jupyter_text>Draw 'fake' data from model: a coin parameter mu and N observed flips<jupyter_code># Define a random generator object, with specific seed for reproducibility
prng = np.random.RandomState(101)

# Settings
N = 10000   # num coin tosses to observe
a = 1      # prior # heads
b = 1      # prior # tails

# Sample a "true" coin frequency from our Beta prior
true_mu = scipy.stats.beta(a,b).rvs(random_state=prng)

# Sample N examples from Bernoulli given mu
x_N = scipy.stats.bernoulli(true_mu).rvs(size=N, random_state=prng)
# OR...
# x_N = prng.rand(N) < true_mu # another equivalent way to draw N samples from a Bernoulli with parameter 'true_mu'<jupyter_output><empty_output><jupyter_text>Exercise 2b: Visualizing Posterior given the *true*, uniform priorComplete the code below, to visualize the posterior over mu for our "fake data" `x_N`<jupyter_code>Npartial_list = [0, 3, 10, 33, 100, 333, 1000, 10000]
L = len(Npartial_list)

fig, axgrid = plt.subplots(
    nrows=1, ncols=L,
    sharex=True, sharey=False,
    figsize=(3 * L, 3))

# Define a grid of mu values at which we'll evaluate the PDF
G = 201
mu_grid_G = np.linspace(0.0001, 0.9999, G)
for aa, Npartial in enumerate(Npartial_list):
    # Compute the parameters of the Beta posterior distribiution, using only first Npartial observations
    ahat, bhat = calc_posterior_parameters(x_N[:Npartial], a, b)

    # Call the PDF function built into scipy.stats.beta, using ahat and bhat as parameters
    # Make sure you evaluate the PDF for each of the G values in the mu_grid_G array
    pdf_grid_G = np.zeros(G) # TODO fix me.
    
    axgrid[aa].plot(mu_grid_G, pdf_grid_G, 'b-')
    axgrid[aa].set_title("after N = %d coins" % Npartial)
plt.tight_layout();<jupyter_output><empty_output><jupyter_text>Discussion 2c : As N gets larger, what happens to the posterior? How does this compare to what happens to the ML estimator? How does this relate to the "true" mu value that generated the data? TODO discuss with your group!  Exercise 2d: Visualizing Posterior given a misinformed priorRepeat Exercise 2b using the *same* fake data, but act as if your prior was instead $a=20, b=3$ (in other words, quite confident in prior belief that the coin would be heads rather than tails.Visualize the posterior over mu for our "fake data" `x_N`, after seeing each of the Npartial examples as in 2b.<jupyter_code>abad = 20
bbad =  3

# TODO adapt code from Exercise 2b, produce a row plot of 7 subplots, one for each Npartial value<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11423076923076923, 'fraction_numerical': 0.04153846153846154, 'mean_word_length': 4.179946879150067, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 30, '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': '242086', 'license': 'MIT', 'n_tokens_mistral': 2916, 'n_tokens_neox': 2617, 'n_words': 912, 'path': 'notebooks/day03-BetaDistribution.ipynb', 'repo_name': 'alex-seto/comp136-21s-assignments', 'size': '7800'}