<jupyter_start><jupyter_text>Chapter 4. Linear Models<jupyter_code>import pandas as pd
import seaborn as sns
import torch

import pyro
import pyro.distributions as dist
import pyro.ops.stats as stats

from rethinking import MAP, coef, extract_samples, link, precis, sim, vcov<jupyter_output><empty_output><jupyter_text>Code 4.1<jupyter_code>pos = torch.empty(1000, 16).uniform_(-1, 1).sum(1)<jupyter_output><empty_output><jupyter_text>Code 4.2<jupyter_code>(1 + torch.empty(12).uniform_(0, 0.1)).prod()<jupyter_output><empty_output><jupyter_text>Code 4.3<jupyter_code>growth = (1 + torch.empty(10000, 12).uniform_(0, 0.1)).prod(1)
sns.distplot(growth, hist=False)
ax = sns.lineplot(growth, dist.Normal(growth.mean(),
                                      growth.std()).log_prob(growth).exp())
ax.lines[1].set_linestyle("--")<jupyter_output><empty_output><jupyter_text>Code 4.4<jupyter_code>big = (1 + torch.empty(10000, 12).uniform_(0, 0.5)).prod(1)
small = (1 + torch.empty(10000, 12).uniform_(0, 0.01)).prod(1)<jupyter_output><empty_output><jupyter_text>Code 4.5<jupyter_code>log_big = (1 + torch.empty(10000, 12).uniform_(0, 0.5)).prod(1).log()<jupyter_output><empty_output><jupyter_text>Code 4.6<jupyter_code>w, n = 6., 9
p_grid = torch.linspace(start=0, end=1, steps=1000)
posterior = (dist.Binomial(n, p_grid).log_prob(torch.tensor(w)).exp()
             * dist.Uniform(0, 1).log_prob(p_grid).exp())
posterior = posterior / posterior.sum()<jupyter_output><empty_output><jupyter_text>Code 4.7<jupyter_code>howell1 = pd.read_csv("../data/Howell1.csv", sep=";")
d = howell1<jupyter_output><empty_output><jupyter_text>Code 4.8<jupyter_code>d.info()
d.head()<jupyter_output><class 'pandas.core.frame.DataFrame'>
RangeIndex: 544 entries, 0 to 543
Data columns (total 4 columns):
height    544 non-null float64
weight    544 non-null float64
age       544 non-null float64
male      544 non-null int64
dtypes: float64(3), int64(1)
memory usage: 17.1 KB<jupyter_text>Code 4.9<jupyter_code>d["height"].head()<jupyter_output><empty_output><jupyter_text>Code 4.10<jupyter_code>d2 = d[d["age"] >= 18]
d2_height = torch.tensor(d2["height"], dtype=torch.float)<jupyter_output><empty_output><jupyter_text>Code 4.11<jupyter_code>x = torch.linspace(100, 250, 101)
sns.lineplot(x, dist.Normal(178, 20).log_prob(x).exp());<jupyter_output><empty_output><jupyter_text>Code 4.12<jupyter_code>x = torch.linspace(-10, 60, 101)
sns.lineplot(x, dist.Uniform(0, 50, validate_args=False).log_prob(x).exp());<jupyter_output><empty_output><jupyter_text>Code 4.13<jupyter_code>sample_mu = torch.empty(int(1e4)).normal_(178, 20)
sample_sigma = torch.empty(int(1e4)).uniform_(0, 50)
prior_h = dist.Normal(sample_mu, sample_sigma).sample()
sns.distplot(prior_h);<jupyter_output><empty_output><jupyter_text>Code 4.14<jupyter_code>mu_list = torch.linspace(start=140, end=160, steps=200)
sigma_list = torch.linspace(start=4, end=9, steps=200)
post = {"mu": mu_list.expand(200, 200).reshape(-1),
        "sigma": sigma_list.expand(200, 200).t().reshape(-1)}
post_LL = dist.Normal(post["mu"],
                      post["sigma"]).log_prob(d2_height.unsqueeze(1)).sum(0)
post_prod = (post_LL + dist.Normal(178, 20).log_prob(post["mu"])
             + dist.Uniform(0, 50).log_prob(post["sigma"]))
post_prob = (post_prod - max(post_prod)).exp()<jupyter_output><empty_output><jupyter_text>Code 4.15<jupyter_code>_, ax = sns.mpl.pyplot.subplots()
ax.contour(post["mu"].reshape(200, 200), post["sigma"].reshape(200, 200),
           post_prob.reshape(200, 200));<jupyter_output><empty_output><jupyter_text>Code 4.16<jupyter_code>_, ax = sns.mpl.pyplot.subplots()
ax.imshow(post_prob.reshape(200, 200), origin="lower",
          extent=(140, 160, 4, 9), aspect="auto")
ax.grid(False)<jupyter_output><empty_output><jupyter_text>Code 4.17<jupyter_code>sample_rows = torch.multinomial(input=post_prob, num_samples=int(1e4),
                                replacement=True)
sample_mu = post["mu"][sample_rows]
sample_sigma = post["sigma"][sample_rows]<jupyter_output><empty_output><jupyter_text>Code 4.18<jupyter_code>ax = sns.scatterplot(sample_mu, sample_sigma, s=64, alpha=0.1, edgecolor="none")
ax.set(xlabel="sample.mu", ylabel="sample.sigma");<jupyter_output><empty_output><jupyter_text>Code 4.19<jupyter_code>sns.distplot(sample_mu)
sns.mpl.pyplot.show()
sns.distplot(sample_sigma);<jupyter_output><empty_output><jupyter_text>Code 4.20<jupyter_code>print(stats.hpdi(sample_mu, 0.89))
print(stats.hpdi(sample_sigma, 0.89))<jupyter_output>tensor([153.8694, 155.1759])
tensor([7.3166, 8.2462])<jupyter_text>Code 4.21<jupyter_code>d3 = stats.resample(d2_height, num_samples=20)<jupyter_output><empty_output><jupyter_text>Code 4.22<jupyter_code>mu_list = torch.linspace(start=150, end=170, steps=200)
sigma_list = torch.linspace(start=4, end=20, steps=200)
post2 = {"mu": mu_list.expand(200, 200).reshape(-1),
         "sigma": sigma_list.expand(200, 200).t().reshape(-1)}
post2_LL = dist.Normal(post2["mu"], post2["sigma"]).log_prob(d3.unsqueeze(1)).sum(0)
post2_prod = (post2_LL + dist.Normal(178, 20).log_prob(post2["mu"])
              + dist.Uniform(0, 50).log_prob(post2["sigma"]))
post2_prob = (post2_prod - max(post2_prod)).exp()
sample2_rows = torch.multinomial(input=post2_prob, num_samples=int(1e4),
                                 replacement=True)
sample2_mu = post2["mu"][sample2_rows]
sample2_sigma = post2["sigma"][sample2_rows]
ax = sns.scatterplot(sample2_mu, sample2_sigma, s=80, alpha=0.1, edgecolor="none")
ax.set(xlabel="mu", ylabel="sigma");<jupyter_output><empty_output><jupyter_text>Code 4.23<jupyter_code>sns.distplot(sample2_sigma, hist=False)
ax = sns.lineplot(sample2_sigma, dist.Normal(sample2_sigma.mean(), sample2_sigma.std())
                  .log_prob(sample2_sigma).exp())
ax.lines[1].set_linestyle("--")<jupyter_output><empty_output><jupyter_text>Code 4.24<jupyter_code>howell1 = pd.read_csv("../data/Howell1.csv", sep=";")
d = howell1
d2 = d[d["age"] >= 18]<jupyter_output><empty_output><jupyter_text>Code 4.25<jupyter_code>def flist(height):
    mu = pyro.sample("mu", dist.Normal(178, 20))
    sigma = pyro.sample("sigma", dist.Uniform(0, 50))
    with pyro.plate("plate"):
        pyro.sample("height", dist.Normal(mu, sigma), obs=height)<jupyter_output><empty_output><jupyter_text>Code 4.26<jupyter_code>d2_height = torch.tensor(d2["height"], dtype=torch.float)
m4_1 = MAP(flist).run(d2_height)<jupyter_output><empty_output><jupyter_text>Code 4.27<jupyter_code>precis(m4_1)<jupyter_output><empty_output><jupyter_text>Code 4.28<jupyter_code>start = {"mu": d2_height.mean(), "sigma": d2_height.std()}<jupyter_output><empty_output><jupyter_text>Code 4.29<jupyter_code>def model(height):
    mu = pyro.sample("mu", dist.Normal(178, 0.1))
    sigma = pyro.sample("sigma", dist.Uniform(0, 50))
    with pyro.plate("plate"):
        pyro.sample("height", dist.Normal(mu, sigma), obs=height)

m4_2 = MAP(model).run(d2_height)
precis(m4_2)<jupyter_output><empty_output><jupyter_text>Code 4.30<jupyter_code>vcov(m4_1)<jupyter_output><empty_output><jupyter_text>Code 4.31<jupyter_code>print(vcov(m4_1).diag())
cov = vcov(m4_1)
print(cov / cov.diag().ger(cov.diag()).sqrt())<jupyter_output>tensor([0.1730, 0.0854])
tensor([[ 1.0000, -0.0025],
        [-0.0025,  1.0000]])<jupyter_text>Code 4.32<jupyter_code>post = extract_samples(m4_1)
{latent: post[latent][:5] for latent in post}<jupyter_output><empty_output><jupyter_text>Code 4.33<jupyter_code>precis(post)<jupyter_output><empty_output><jupyter_text>Code 4.34<jupyter_code>post = dist.MultivariateNormal(torch.stack(list(coef(m4_1).values())),
                               vcov(m4_1)).sample(torch.Size([int(1e4)]))<jupyter_output><empty_output><jupyter_text>Code 4.35<jupyter_code>def model(height):
    mu = pyro.sample("mu", dist.Normal(178, 20))
    log_sigma = pyro.sample("log_sigma", dist.Normal(2, 10))
    with pyro.plate("plate"):
        pyro.sample("height", dist.Normal(mu, log_sigma.exp()), obs=height)

m4_1_logsigma = MAP(model).run(d2_height)<jupyter_output><empty_output><jupyter_text>Code 4.36<jupyter_code>post = extract_samples(m4_1_logsigma)
sigma = post["log_sigma"].exp()<jupyter_output><empty_output><jupyter_text>Code 4.37<jupyter_code>sns.scatterplot("weight", "height", data=d2);<jupyter_output><empty_output><jupyter_text>Code 4.38<jupyter_code># load data again, since it's a long way back
howell1 = pd.read_csv("../data/Howell1.csv", sep=";")
d = howell1
d2 = d[d["age"] >= 18]

# fit model
def model(weight, height):
    a = pyro.sample("a", dist.Normal(178, 100))
    b = pyro.sample("b", dist.Normal(0, 10))
    mu = a + b * weight
    sigma = pyro.sample("sigma", dist.Uniform(0, 50))
    with pyro.plate("plate"):
        pyro.sample("height", dist.Normal(mu, sigma), obs=height)

d2_weight = torch.tensor(d2["weight"], dtype=torch.float)
d2_height = torch.tensor(d2["height"], dtype=torch.float)
m4_3 = MAP(model).run(d2_weight, d2_height)<jupyter_output><empty_output><jupyter_text>Code 4.39<jupyter_code>def model(weight, height):
    a = pyro.sample("a", dist.Normal(178, 100))
    b = pyro.sample("b", dist.Normal(0, 10))
    sigma = pyro.sample("sigma", dist.Uniform(0, 50))
    with pyro.plate("plate"):
        pyro.sample("height", dist.Normal(a + b * weight, sigma), obs=height)

m4_3 = MAP(model).run(d2_weight, d2_height)<jupyter_output><empty_output><jupyter_text>Code 4.40<jupyter_code>precis(m4_3)<jupyter_output><empty_output><jupyter_text>Code 4.41<jupyter_code>precis(m4_3, corr=True)<jupyter_output><empty_output><jupyter_text>Code 4.42<jupyter_code>d2_weight_c = d2_weight - d2_weight.mean()<jupyter_output><empty_output><jupyter_text>Code 4.43<jupyter_code>m4_4 = MAP(model).run(d2_weight_c, d2_height)<jupyter_output><empty_output><jupyter_text>Code 4.44<jupyter_code>precis(m4_4, corr=True)<jupyter_output><empty_output><jupyter_text>Code 4.45<jupyter_code>sns.scatterplot("weight", "height", data=d2)
x = torch.linspace(30, 65, 101)
sns.lineplot(x, (coef(m4_3)["a"] + coef(m4_3)["b"] * x), color="k");<jupyter_output><empty_output><jupyter_text>Code 4.46<jupyter_code>post = extract_samples(m4_3)<jupyter_output><empty_output><jupyter_text>Code 4.47<jupyter_code>{latent: post[latent][:5].detach() for latent in post}<jupyter_output><empty_output><jupyter_text>Code 4.48<jupyter_code>N = 10
dN = {"weight": d2_weight[:N], "height": d2_height[:N]}
mN = MAP(model).run(**dN)<jupyter_output><empty_output><jupyter_text>Code 4.49<jupyter_code># extract 20 samples from the posterior
idx = mN._categorical.sample(torch.Size([20]))
post = {latent: samples[idx] for latent, samples in extract_samples(mN).items()}

# display raw data and sample size
ax = sns.scatterplot("weight", "height", data=dN)
ax.set(xlabel="weight", ylabel="height", title="N = {}".format(N))

# plot the lines, with transparency
x = torch.linspace(30, 65, 101)
for i in range(20):
    sns.lineplot(x, post["a"][i] + post["b"][i] * x, color="k", alpha=0.3)<jupyter_output><empty_output><jupyter_text>Code 4.50<jupyter_code>post = extract_samples(m4_3)
mu_at_50 = post["a"] + post["b"] * 50<jupyter_output><empty_output><jupyter_text>Code 4.51<jupyter_code>ax = sns.distplot(mu_at_50)
ax.set(xlabel="mu|weight=50", ylabel="Density");<jupyter_output><empty_output><jupyter_text>Code 4.52<jupyter_code>stats.hpdi(mu_at_50, prob=0.89)<jupyter_output><empty_output><jupyter_text>Code 4.53<jupyter_code>mu = link(m4_3)
mu.shape, mu[:5, 0]<jupyter_output><empty_output><jupyter_text>Code 4.54<jupyter_code># define sequence of weights to compute predictions for
# these values will be on the horizontal axis
weight_seq = torch.arange(start=25., end=71, step=1)

# use link to compute mu
# for each sample from posterior
# and for each weight in weight_seq
mu = link(m4_3, data={"weight": weight_seq})
mu.shape, mu[:5, 0]<jupyter_output><empty_output><jupyter_text>Code 4.55<jupyter_code># use visible=False to hide raw data
sns.scatterplot("weight", "height", data=d2, visible=False)

# loop over samples and plot each mu value
for i in range(100):
    sns.scatterplot(weight_seq, mu[i], color="royalblue", alpha=0.1)<jupyter_output><empty_output><jupyter_text>Code 4.56<jupyter_code># summarize the distribution of mu
mu_mean = mu.mean(0)
mu_HPDI = stats.hpdi(mu, prob=0.89, dim=0)<jupyter_output><empty_output><jupyter_text>Code 4.57<jupyter_code># plot raw data
# fading out points to make line and interval more visible
sns.scatterplot("weight", "height", data=d2, alpha=0.5)

# plot the MAP line, aka the mean mu for each weight
ax = sns.lineplot(weight_seq, mu_mean, color="k")

# plot a shaded region for 89% HPDI
ax.fill_between(weight_seq, mu_HPDI[0], mu_HPDI[1], color="k", alpha=0.2);<jupyter_output><empty_output><jupyter_text>Code 4.58<jupyter_code>post = extract_samples(m4_3)
mu_link = lambda weight: post["a"].unsqueeze(1) + post["b"].unsqueeze(1) * weight
weight_seq = torch.arange(start=25., end=71, step=1)
mu = mu_link(weight_seq)
mu_mean = mu.mean(0)
mu_HPDI = stats.hpdi(mu, prob=0.89, dim=0)<jupyter_output><empty_output><jupyter_text>Code 4.59<jupyter_code>sim_height = sim(m4_3, data={"weight": weight_seq})
sim_height.shape, sim_height[:5, 0]<jupyter_output><empty_output><jupyter_text>Code 4.60<jupyter_code>height_PI = stats.pi(sim_height, prob=0.89, dim=0)<jupyter_output><empty_output><jupyter_text>Code 4.61<jupyter_code># plot raw data
sns.scatterplot("weight", "height", data=d2, alpha=0.5)

# draw MAP line
ax = sns.lineplot(weight_seq, mu_mean, color="k")

# draw HPDI region for line
ax.fill_between(weight_seq, mu_HPDI[0], mu_HPDI[1], color="k", alpha=0.15)

# draw PI region for simulated heights
ax.fill_between(weight_seq, height_PI[0], height_PI[1], color="k", alpha=0.15);<jupyter_output><empty_output><jupyter_text>Code 4.62<jupyter_code>sim_height = sim(m4_3, data={"weight": weight_seq}, n=int(1e4))
height_PI = stats.pi(sim_height, prob=0.89, dim=0)<jupyter_output><empty_output><jupyter_text>Code 4.63<jupyter_code>def sim_fn(weight):
    mean = post["a"].unsqueeze(1) + post["b"].unsqueeze(1) * weight
    sd = post["sigma"].unsqueeze(1)
    return dist.Normal(loc=mean, scale=sd).sample()

post = extract_samples(m4_3)
weight_seq = torch.arange(start=25., end=71, step=1)
sim_height = sim_fn(weight_seq)
height_PI = stats.pi(sim_height, prob=0.89, dim=0)<jupyter_output><empty_output><jupyter_text>Code 4.64<jupyter_code>howell1 = pd.read_csv("../data/Howell1.csv", sep=";")
d = howell1
d.info()
d.head()<jupyter_output><class 'pandas.core.frame.DataFrame'>
RangeIndex: 544 entries, 0 to 543
Data columns (total 4 columns):
height    544 non-null float64
weight    544 non-null float64
age       544 non-null float64
male      544 non-null int64
dtypes: float64(3), int64(1)
memory usage: 17.1 KB<jupyter_text>Code 4.65<jupyter_code>weight = torch.tensor(d["weight"], dtype=torch.float)
weight_s = (weight - weight.mean()) / weight.std()<jupyter_output><empty_output><jupyter_text>Code 4.66<jupyter_code>def model(weight, weight2, height):
    a = pyro.sample("a", dist.Normal(178, 100))
    b1 = pyro.sample("b1", dist.Normal(0, 10))
    b2 = pyro.sample("b2", dist.Normal(0, 10))
    mu = a + b1 * weight + b2 * weight2
    sigma = pyro.sample("sigma", dist.Uniform(0, 50))
    with pyro.plate("plate"):
        pyro.sample("height", dist.Normal(mu, sigma), obs=height)

weight_s2 = weight_s ** 2
height = torch.tensor(d["height"], dtype=torch.float)
m4_5 = MAP(model).run(weight_s, weight_s2, height)<jupyter_output><empty_output><jupyter_text>Code 4.67<jupyter_code>precis(m4_5)<jupyter_output><empty_output><jupyter_text>Code 4.68<jupyter_code>weight_seq = torch.linspace(start=-2.2, end=2, steps=30)
pred_data = {"weight": weight_seq, "weight2": weight_seq ** 2}
mu = link(m4_5, data=pred_data)
mu_mean = mu.mean(0)
mu_PI = stats.pi(mu, prob=0.89, dim=0)
sim_height = sim(m4_5, data=pred_data)
height_PI = stats.pi(sim_height, prob=0.89, dim=0)<jupyter_output><empty_output><jupyter_text>Code 4.69<jupyter_code>ax = sns.scatterplot(weight_s, height, alpha=0.5)
ax.set(xlabel="weight.s", ylabel="height")
sns.lineplot(weight_seq, mu_mean, color="k")
ax.fill_between(weight_seq, mu_PI[0], mu_PI[1], color="k", alpha=0.2)
ax.fill_between(weight_seq, height_PI[0], height_PI[1], color="k", alpha=0.2);<jupyter_output><empty_output><jupyter_text>Code 4.70<jupyter_code>def model(weight, weight2, weight3, height):
    a = pyro.sample("a", dist.Normal(178, 100))
    b1 = pyro.sample("b1", dist.Normal(0, 10))
    b2 = pyro.sample("b2", dist.Normal(0, 10))
    b3 = pyro.sample("b3", dist.Normal(0, 10))
    mu = a + b1 * weight + b2 * weight2 + b3 * weight3
    sigma = pyro.sample("sigma", dist.Uniform(0, 50))
    with pyro.plate("plate"):
        pyro.sample("height", dist.Normal(mu, sigma), obs=height)

weight_s3 = weight_s ** 3
m4_6 = MAP(model).run(weight_s, weight_s2, weight_s3, height)<jupyter_output><empty_output><jupyter_text>Code 4.71<jupyter_code>fig, ax = sns.mpl.pyplot.subplots()
sns.scatterplot(weight_s, height, alpha=0.5)
ax.set(xlabel="weight", ylabel="height", xticks=[]);<jupyter_output><empty_output><jupyter_text>Code 4.72<jupyter_code>at = torch.tensor([-2, -1, 0, 1, 2])
labels = at * weight.std() + weight.mean()
ax.set_xticks(at)
ax.set_xticklabels([round(label.item(), 1) for label in labels])
fig<jupyter_output><empty_output><jupyter_text>Code 4.73<jupyter_code>sns.scatterplot("weight", "height", data=howell1, alpha=0.4);<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup