How do I find where a Dimensionality error is when I can't graph my model?

[jordan.howell2]

Hello,

I’m getting a dimensionality error trying to graph my model but can’t find where the error is.

---------------------------------------------------------------------------
ShapeError                                Traceback (most recent call last)
/tmp/ipykernel_11866/2705921953.py in <module>
     86               predictions,
     87               observed=eaches,
---> 88                dims='obs_id'
     89         )
     90 

/opt/conda/lib/python3.7/site-packages/pymc/distributions/distribution.py in __new__(cls, name, *args, **kwargs)
    592             raise ValueError("Transformations for discrete distributions")
    593 
--> 594         return super().__new__(cls, name, *args, **kwargs)
    595 
    596 

/opt/conda/lib/python3.7/site-packages/pymc/distributions/distribution.py in __new__(cls, name, rng, dims, initval, observed, total_size, transform, *args, **kwargs)
    276             dims=dims,
    277             transform=transform,
--> 278             initval=initval,
    279         )
    280 

/opt/conda/lib/python3.7/site-packages/pymc/model.py in register_rv(self, rv_var, name, data, total_size, dims, transform, initval)
   1311             # `rv_var` is potentially changed by `make_obs_var`,
   1312             # for example into a new graph for imputation of missing data.
-> 1313             rv_var = self.make_obs_var(rv_var, data, dims, transform)
   1314 
   1315         return rv_var

/opt/conda/lib/python3.7/site-packages/pymc/model.py in make_obs_var(self, rv_var, data, dims, transform)
   1338         if data.ndim != rv_var.ndim:
   1339             raise ShapeError(
-> 1340                 "Dimensionality of data and RV don't match.", actual=data.ndim, expected=rv_var.ndim
   1341             )
   1342 

ShapeError: Dimensionality of data and RV don't match. (actual 1 != expected 2)

When I try to insert a delta.shape.eval() or m.shape.eval() I then get the following error.

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
/tmp/ipykernel_11866/4216889251.py in <module>
     52     mu_delta = pm.Laplace('mu_delta', 0,1, shape=(n_groups, n_changepoints), dims='location')
     53     delta = pm.Deterministic('delta', sigma_delta*mu_delta)
---> 54     delta.shape_eval()
     55     growth = ((k[group] + at.dot(a_, delta[group])))
     56     offset = ((m[group] + at.dot(a_, -s_ * delta[group])))

AttributeError: 'TensorVariable' object has no attribute 'shape_eval'

For reference, here is my model code.

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

group = df_train['location'].cat.codes.values
group_mapping = dict(enumerate(df_train['location'].cat.categories))
n_groups = df_train['location'].nunique()

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)
    eaches = pm.MutableData('eaches', df_train['eaches'])
    
    mu_k = pm.Normal('mu_k', 0,100, shape = n_groups, dims = 'location')
    sigma_k = pm.HalfNormal('sigma_k', 5)
    k = pm.Deterministic('k', mu_k*sigma_k)

    m = pm.Normal('m', 0, 5, shape = n_groups)

    sigma_delta = pm.HalfCauchy('sigma_delta', .05)
    mu_delta = pm.Laplace('mu_delta', 0,1, shape=(n_groups, n_changepoints), dims='location')
    delta = pm.Deterministic('delta', sigma_delta*mu_delta)

    growth = ((k[group] + at.dot(a_, delta[group])))
    offset = ((m[group] + at.dot(a_, -s_ * delta[group])))
    trend = pm.Deterministic('trend', growth * t_ + offset)
    
    predictions =  pm.Deterministic('predictions', np.exp(trend[location_idx])) 

    pm.Poisson('obs',
              predictions,
              observed=eaches,
               dims='obs_id'
        )

    
pm.model_to_graphviz( partial_pooled_model)

Code adapted from @MBrouns repo timeseers/timeseries_model.py at master · MBrouns/timeseers (github.com)

UPDATE:
I took out the PM.Poisson... code at the end and got it to graph this:

I’m not sure how trend is ending up with 715x715 shape or how to fix it.

[jordan.howell2]

I figured it out. I had to replace the growth/offset code with pm.math instead of at.dot.

New Model Code:

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,100, shape = n_groups, dims = 'location')
    sigma_k = pm.HalfCauchy('sigma_k', 1)
    k = pm.Deterministic('k', mu_k*sigma_k)

    m = pm.Normal('m', 0, 5, shape=n_groups)

    sigma_delta = pm.HalfCauchy('sigma_delta', .05)
    mu_delta = pm.Laplace('mu_delta', 0,1, shape=(n_groups, n_changepoints), dims='location')
    delta = pm.Deterministic('delta', sigma_delta*mu_delta)

    growth = pm.Deterministic('growth',  (k[group] + pm.math.sum(a_ * delta[group], axis=1)))
    gamma = -s_ * delta[group, :]
    offset = pm.Deterministic('offset', (m[group] + pm.math.sum(a_ * gamma, axis=1)))
    trend = pm.Deterministic('trend', growth * t_ + offset)

    pm.Poisson('obs',
                  predictions,
                  observed=eaches,
                   dims='obs_id'
            

I’ll have to admit, I don’t understand “why” I have to do this. It works, but I don’t really feel like I’m learning in this instance.

[drbenvincent]

I’ll have to admit, I don’t understand “why” I have to do this. It works, but I don’t really feel like I’m learning in this instance.

Hi! Have a look at the graphs side by side. Forgetting the non-important changes in the visual aspects, these are not identical models. In the first one trend is shape 715x715, in the second one it is 715.

So these are doing fundamentally different things, and I suspect the first model was not doing what you wanted it do.

In terms of learning, I’d also suggest to really embrace the use of dims. I see that you’ve got partial use of dims at the moment… some of the plates in the graph digram have dimension names, but others don’t. My suggestion would be to use dims more (and move away from shape) and to really inspect the graphviz for correctness It can take a while to spot potential errors, but it is useful.

Hope that helps in terms of your frustration of not learning from that particular solution you found.

[jordan.howell2]

Thank you for taking the time on a thoughtful answer.

So if I use dims, I don’t need to use shape?

[cluhmann]

@jordan.howell2

So if I use dims, I don’t need to use shape?

Correct. You can check out Oriol’s post and my own, both of which touch on the use of coords to provide shape info in a more user-friendly manner.

[jordan.howell2]

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.

[cluhmann]

@jordan.howell2

But what about instances where a parameter has more than on dimension in the shape?

This is fine. Coords and dims work in multiple dimensions:

coords = {'day':np.arange(7),
          'city':['Berlin', 'London', 'Tokyo']
         }
with pm.Model(coords=coords) as model:
    mu = pm.Normal('mu', 0, 1)
    # shape
    #normals = pm.Normal('normals', mu, 1, shape=(7, 3))
    # dims
    normals = pm.Normal('normals', mu, 1, dims=('day', 'city'))

    idata = pm.sample()

But I don’t see any instances of multidimensional parameter arrays in your model (though I may have missed it).

As for the shape error you’re getting, it’s not clear to me what is happening because many things have changed between the two models.

[jordan.howell2]

‘’’
As for the shape error you’re getting, it’s not clear to me what is happening because many things have changed between the two models.
‘’’

Unfortunately that part is true. It’s a result of me trying anything to get this to work. I’m starting over with the non-hierarchical model and making one piece at a time hierarchical. I think I have trend working for the hierarchy. I think at this point, the question is good. If something comes up at this point, with the model so different, I’ll repost under a different thread. Thank you for the posts. Learning slowly but surely!