In the GLM: Hierarchical Linear Regression article a hierarchical regression model is used to predict radon levels in different counties.
pm.sample the parameter’s distributions for each country can be fitted. Now that I know the parameters, how can I sample the model for a specific county only?
Full details. Code from article:
import matplotlib.pyplot as plt import numpy as np import pymc3 as pm import pandas as pd import theano data = pd.read_csv(pm.get_data('radon.csv')) data['log_radon'] = data['log_radon'].astype(theano.config.floatX) county_names = data.county.unique() county_idx = data.county_code.values n_counties = len(data.county.unique()) with pm.Model() as hierarchical_model: # Hyperpriors for group nodes mu_a = pm.Normal('mu_a', mu=0., sd=100**2) sigma_a = pm.HalfCauchy('sigma_a', 5) mu_b = pm.Normal('mu_b', mu=0., sd=100**2) sigma_b = pm.HalfCauchy('sigma_b', 5) # Intercept for each county, distributed around group mean mu_a # Above we just set mu and sd to a fixed value while here we # plug in a common group distribution for all a and b (which are # vectors of length n_counties). a = pm.Normal('a', mu=mu_a, sd=sigma_a, shape=n_counties) # Intercept for each county, distributed around group mean mu_a b = pm.Normal('b', mu=mu_b, sd=sigma_b, shape=n_counties) # Model error eps = pm.HalfCauchy('eps', 5) radon_est = a[county_idx] + b[county_idx] * data.floor.values # Data likelihood radon_like = pm.Normal('radon_like', mu=radon_est, sd=eps, observed=data.log_radon)
Then to fit the model and determine parameter distributions:
with hierarchical_model: hierarchical_trace = pm.sample(draws=2000, n_init=1000)
How can I now sample the model for a specific county only, for example AITKIN?