Radon Hierarchical Model with state as a level in addition to county. Causes error in PyMC3 of "fatal error: bracket nesting level exceeded maximum of 256"

foo.py

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pymc3 as pm 
import pandas as pd

data = pd.read_csv('./radon.csv')

state_names = data.state.unique()
county_names = data.county.unique()
state_county_names = set(list(zip(data.state.values, data.county.values)))

with pm.Model() as hierarchical_model:
    global_m = pm.Normal('global_m', mu=0., sd=100**2)
    global_m_sigma = pm.Uniform('global_m_sigma', lower=0, upper=100)
    global_b = pm.Normal('global_b', mu=0., sd=100**2)
    global_b_sigma = pm.Uniform('global_b_sigma', lower=0, upper=100)
    
    states = {state_name: {'mu': pm.Normal(state_name + ' mu', mu=global_m, sd=global_m_sigma), 'b': pm.Normal(state_name + ' b', mu=global_b, sd=global_b_sigma)} for state_name in state_names}
    counties = {}
    for state_name, county_name in state_county_names:
      county_df = data[(data['state'] == state_name) & (data['county'] == county_name)]
      county_slope = pm.Normal('{0}_{1}_mu'.format(state_name, county_name), mu=states[state_name]['mu'], sd=pm.Uniform('{0}_{1}_m_sigma'.format(state_name, county_name), lower=0, upper=100))
      county_intercept = pm.Normal('{0}_{1}_b'.format(state_name, county_name), mu=states[state_name]['b'], sd=pm.Uniform('{0}_{1}_b_sigma'.format(state_name, county_name), lower=0, upper=100))
      county_prediction = county_slope * county_df['floor'].values + county_intercept
      county_prediction_error = pm.Normal('{0}_{1}_error'.format(state_name, county_name), mu=0.0, sd=county_prediction - county_df['log_radon'].values)
      counties[(state_name, county_name)] = {'mu': county_slope, 'b': county_intercept, 'prediction': county_prediction_error}

with hierarchical_model:
    hierarchical_trace = pm.sample(10000)