Bayesian regression confidence intervals with Pymc3 This question is based on question 1 of the week 2 Statistical Rethinking problems, i.e. q1 here:  https://github.com/rmcelreath/stat_rethinking_2022/blob/main/homework/week02.pdf
I have a pandas data frame df with observed weights and heights (df.weight and df.height resp.)
I am trying to use linear regression (predicting weight from height) with Pymc3 to find an 89% confidence interval for the weight of an 140cm tall individual.
This is my setup so far
with pm.Model() as q1_model:
    # specify the model 
    alpha = pm.Normal('alpha', mu=45, sd=100)
    beta = pm.Normal('beta', mu=0, sd=10)
    sigma = pm.Uniform('sigma', lower=0, upper=50)
    weight = pm.Normal('weight', mu=alpha + beta * df.height_c, sd=sigma, observed=df.weight)
    
    # find the posterior distribution of the weight
    trace = pm.sample(1000)

where df.height_c is df.height - df.height.mean(). But I'm not sure how to get the interval. I tried:
pm.hdi(trace.alpha + (140 - df.height.mean()) * trace.beta, 0.89)

which gives [35.0, 36.6] but I think this is just an interval for the mean height mu of an 140cm person, not for height itself. The actual interval is supposedly meant to be [29.1, 42.8] confirming my suspicion here.

This exercise asks for a model of weight given height and a 89% compatibility interval for the weight of an individual with height = 140cm.
Aside: @GrahamBornholt comments that "it is a credible interval you wish to calculate". So let's clarify the terminology. McElreath prefers the term "compatibility interval" to highlight that the "the interval indicates a range of parameter values compatible with the model and data." The crucial thing to keep in mind is not the name of the interval but its interpretation, which is fundamentally different from the meaning of a frequentist confidence interval. See What's the difference between a confidence interval and a credible interval?
I'll use a bit of math to go with the python code. Let $x$ denote height (in cm) and $y$ denote weight (in kg). Then the regression model for $y$ given $x$ is:
$$
\begin{aligned}
y = \alpha + \beta(x - \bar{x}) + e, \quad e \sim \operatorname{N}(0,\sigma^2)
\end{aligned}
$$
The weight of a 140cm-tall individual ($x = 140$) is modeled as the sum of two quantities: the average weight of 140cm-tall persons $\operatorname{E}(y) = \alpha + \beta(140 - \bar{x})$ plus the difference $e$ between the individual's weight and the average weight at that height; this difference is called the error.
So to construct a compatibility interval for the weight of a person (rather than a compatibility interval for the average weight) at a specific height, you have to simulate both $\operatorname{E}(y)$ and the error from the posterior. As McElreath explains in Section 4.4.3.5: sample from the posterior to simulate $y$s, not a distribution of plausible average $y$.
First set up and fit the model.
import numpy as np
import pandas as pd
import pymc3 as pm
import scipy.stats as stats

np.random.seed(1234)

howell = pd.read_csv("Howell1.csv", sep=";")
adults = howell.query("age >= 18")
x = adults["height"].values
y = adults["weight"].values
xbar = np.mean(x)

with pm.Model() as model:
    alpha = pm.Normal("alpha", mu=45, sd=100)
    beta = pm.Normal("beta", mu=0, sigma=10)
    sigma = pm.Uniform("sigma", lower=0, upper=50)
    weight = pm.Normal("weight", mu=alpha + beta * (x - xbar), sd=sigma, observed=y)
    trace = pm.sample(1000, tune=1000)

The exercise asks about the weight of a 140cm-tall person.
x0 = 140

First we construct a 89% compatibility interval to represent the uncertainty about the average weight of 140cm-tall individuals.
pp_Ey_at_x0 = trace.alpha + trace.beta * (x0 - xbar)
pm.hdi(pp_Ey_at_x0, hdi_prod=0.89)
#> [34.8, 36.7]

Then we construct a 89% compatibility interval to represent the uncertainty about the weight of a person who is 140cm tall. [This is also called a prediction interval.]
pp_y_at_x0 = pp_Ey_at_x0 + stats.norm.rvs(0, trace.sigma)
pm.hdi(pp_y_at_x0, hdi_prod=0.89)
#> [27.8, 43.6]

PS. It would be worth it to revisit the prior specification. Does the prior weight ~ $\operatorname{N}(45,100)$ at the average adult height seem plausible? Careful consideration about the prior is a big part of Statistical Rethinking. And it's a big part of doing Bayesian analysis well. See Section 4.4.1.3 on choosing priors and especially Figure 4.5 for an illustration of prior predictive simulation.
