Thank you. I read both. But what about instances where a parameter has more than on dimension in the shape?
For example, here is a non-hierarchical model. Notice the shape on delta and two seasonality parameters (yearly and monthly beta).
n_changepoints = 8
t = np.arange(len(df_sat)) / len(df_sat)
s = np.linspace(0, np.max(t), n_changepoints+2)[1:-1]
a = (t[:, None] > s)*1
n_components = 10
def yearly_X(t, p=365.25, n=n_components):
x = 2 * np.pi * (np.arange(n)+1) * t[:, None] / p
return np.concatenate((np.cos(x), np.sin(x)), axis = 1)
monthly_n_components = 5
def monthly_X(t, p=30.5, n=monthly_n_components):
x = 2 * np.pi * (np.arange(n)+1) * t[:, None] / p
return np.concatenate((np.cos(x), np.sin(x)), axis = 1)
# df_sat['promo_flag_month_before'] = df_sat['promo_flag'].shift(1)
# df_sat['promo_flag_month_before'] = df_sat['promo_flag_month_before'].fillna(method = 'bfill')
# promo_before = np.array(df_sat['promo_flag_month_before'])
# promo = np.array(df_sat['promo_flag'])
# month = np.array(df_sat['month'])
yearly_fourier = yearly_X(t, 365.25/len(t))
monthly_fourier = monthly_X(t, 30.5/len(t))
with pm.Model() as model:
# promo_before = pm.MutableData('promo_before', promo_before)
# promo = pm.MutableData('promo', promo
# month_ = pm.MutableData('month', month)
t_ = pm.MutableData('t', t)
a_ = pm.MutableData('a', a)
s_ = pm.MutableData('s', s)
yearly_f_ = pm.MutableData('yearly_f_', yearly_fourier)
monthly_m_ = pm.MutableData('monthly_m_',monthly_fourier)
sept_month = pm.MutableData('sept_month',sept)
march_month = pm.MutableData('march_month', march)
eaches = pm.MutableData('eaches', df_sat['eaches'])
k = pm.Normal('k', 0, 1)
m = pm.Normal('m', 0, 5)
delta = pm.Laplace('delta', 0, 0.05, shape=n_changepoints)
growth = k + at.dot(a_, delta)
offset = (m + at.dot(a_, -s_ * delta))
trend = pm.Deterministic('trend', growth * t_ + offset)
yearly_beta = pm.Normal('yearly_beta', 0, 1, shape = n_components*2)
yearly_seasonality = pm.Deterministic('yearly_seasonality',at.dot(yearly_f_, yearly_beta))
monthly_beta = pm.Normal('monthly_beta', 0, 10, shape = monthly_n_components*2)
monthly_seasonality = pm.Deterministic('monthly_seasonality',at.dot(monthly_m_, monthly_beta))
sept_beta = pm.Normal('sept_beta', 0,10)
march_beta = pm.Normal('march_beta', 0,10)
# promo_beta_before = pm.Normal('promo_beta_before', 0, 1)
# promo_beta = pm.Normal('promo_beta', 0, 1)
# promo_beta_interaction = pm.Normal('promo_beta_interaction', 0, 1)
predictions = pm.Deterministic('predictions', np.exp(trend + monthly_seasonality + yearly_seasonality + (sept_beta*sept_month) + (march_beta*march_month)))
# error = pm.HalfCauchy('error', 0.5)
pm.Poisson('obs',
predictions,
observed=eaches
)
# trace = pymc.sampling_jax.sample_numpyro_nuts(tune=2000, draws = 2000)
When I try to do this in a hierarchical using dims, I’m getting a dimension error.
n_changepoints = 8
t = np.arange(len(df_train)) / len(df_train)
s = np.linspace(0, np.max(t), n_changepoints+2)[1:-1]
a = (t[:, None] > s)*1
n_components = 10
def yearly_X(t, p=365.25, n=n_components):
x = 2 * np.pi * (np.arange(n)+1) * t[:, None] / p
return np.concatenate((np.cos(x), np.sin(x)), axis = 1)
monthly_n_components = 5
def monthly_X(t, p=30.5, n=monthly_n_components):
x = 2 * np.pi * (np.arange(n)+1) * t[:, None] / p
return np.concatenate((np.cos(x), np.sin(x)), axis = 1)
month = pd.get_dummies(df_train['month'])
sept = np.array(month[9])
march = np.array(month[3])
month = np.array(df_train['month'])
yearly_fourier = yearly_X(t, 365.25/len(t))
monthly_fourier = monthly_X(t, 30.5/len(t))
location_idxs, locations = pd.factorize(df_train['location'])
coords_={
"location":df_train['location'].unique(),
"obs_id":np.arange(len(location_idxs))
}
with pm.Model(coords=coords_) as partial_pooled_model:
location_idx = pm.ConstantData('location_idx', location_idxs, dims = 'obs_id')
t_ = pm.MutableData('t', t)
a_ = pm.MutableData('a', a)
s_ = pm.MutableData('s', s)
# yearly_f_ = pm.MutableData('yearly_f_', yearly_fourier)
# monthly_m_ = pm.MutableData('monthly_m_',monthly_fourier)
# sept_month = pm.MutableData('sept_month',sept)
# march_month = pm.MutableData('march_month', march)
# eaches = pm.MutableData('eaches', df_train['eaches'])
mu_k = pm.Normal('mu_k', 0, 5,)
sigma_k = pm.HalfCauchy('sigma_k', 5)
k = pm.Normal('k', mu_k, sigma_k, dims = 'location')
m = pm.Normal('m', 0, 5, dims = 'location')
mu_delta = pm.Normal('mu_delta', 0, 1, dims = 'location')
sigma_delta = pm.HalfCauchy('sigma_delta', .05)
delta = pm.Laplace('delta', mu_delta, sigma_delta, dims = 'location')
growth = pm.Deterministic('growth', k[location_idx] + at.dot(a_, delta[location_idx]), dims = 'location')
offset = pm.Deterministic('offset', (m[location_idx] + at.dot(a_, (-s_ * delta[location_idx]))), dims = 'location')
trend = pm.Deterministic('trend', growth[location_idx] * t_ + offset[location_idx], dims = 'location')
# mu_yearly = pm.Normal('mu_yearly', 0, 5)
# sigma_yealry = pm.HalfCauchy('sigma_yearly', 5)
# yearly_beta = pm.Normal('yearly_beta', mu_yearly, 1, dims = 'location')
# yearly_seasonality = pm.Deterministic('yearly_seasonality',at.dot(yearly_f_, yearly_beta[location_idx]))
# mu_monthly = pm.Normal('mu_monthly',0,5)
# sigma_monthly = pm.HalfCauchy('sigma_monthly',10)
# monthly_beta = pm.Normal('monthly_beta', mu_monthly, sigma_monthly, dims = 'location')
# monthly_seasonality = pm.Deterministic('monthly_seasonality',at.dot(monthly_m_, monthly_beta[location_idx]))
# sept_beta = pm.Normal('sept_beta', 0,10)
# march_beta = pm.Normal('march_beta', 0,10)
# promo_beta_before = pm.Normal('promo_beta_before', 0, 1),
# promo_beta = pm.Normal('promo_beta', 0, 1)
# promo_beta_interaction = pm.Normal('promo_beta_interaction', 0, 1)
predictions = pm.Deterministic('predictions', np.exp(trend[location_idx]))# + monthly_seasonality[location_idx] + yearly_seasonality[location_idx] + (sept_beta*sept_month) + (march_beta*march_month)))
# error = pm.HalfCauchy('error', 0.5)
pm.Poisson('obs',
predictions,
observed=eaches
)
# trace = pymc.sampling_jax.sample_numpyro_nuts(tune=2000, draws = 2000)
pm.model_to_graphviz( partial_pooled_model
When I try to sample this, I get
ValueError: shapes (715,8) and (715,) not aligned: 8 (dim 1) != 715 (dim 0)
This tells me the problem lies at a…or maybe delta but I’m not sure how to make it work using dims.