I am working on the same dataset as you following along the book and implementing the model using Pymc. Using the data of the divorce I am getting some really bad sampling but in the book it doesn't seem to be a problem? do you have any idea why? I am trying with your code as we speak.
this is my model, my values are standardized
my model:
with pm.Model(coords=coords_div) as Model_missing_data:
a = pm.Normal("intercept",0,0.2)
bA = pm.Normal("beta_age",0,0.5)
bM = pm.Normal("beta_marriage",0,0.5)
mu_div = a + bA*age + bM*mar
sigma_div = pm.Exponential("error",1)
# true divorce rate
div_true = pm.Normal("Divorce_True",mu_div, sigma_div,dims="obs_id")
# Observed divorce rate generated by the (unknown) true div rate
D_obs = pm.Normal("Divorce_Obs",mu = div_true, sigma =div_sd ,observed=div)
```
thanks for your work,