Issue with dimension mismatch, unsure about index variable behavior

I’m trying to figure out why I am having a dimension mismatch here. It was my understanding that we could access the index of a shared variable to pass it along to another process. However, the line where “kappa” is declared kicks off an error due to a dimension mismatch. For some reason m1[build_idx] has dimension 64 rather than a dimension of 1 as expected. What am I missing? I could not find a good example in the documentation nor in any other reported issue. Thanks in advance!

Y = observations['z']
X = observations['x']
build_idx, build = pd.factorize(observations["building"], sort=True)
group_idx,group = pd.factorize(observations['group'],sort=True)
coords = {"group":group,"build":build,"obs_id":range(observations.shape[0])}
with pm.Model(coords=coords) as model:
    
    # Setup indexing variables
    build_idx = pm.Data("build_idx", build_idx, dims="obs_id")
    group_idx = pm.Data("group_idx", group_idx, dims="obs_id")
    x = pm.Data("x", X)
    
    # Hyper-priors
    m1 = pm.Uniform('m1',0,10, dims = 'build')
    m2 = pm.Uniform('m2',0,10, dims = 'build')
    s1 = pm.Uniform('s1',0,2, dims = 'build')
    s2 = pm.Uniform('s2',0,2, dims = 'build')
        
    kappa = pm.Normal("kappa", m1[build_idx], s1[build_idx], dims = "group")
    llambda = pm.Normal("llambda",m2[build_idx], s2[build_idx], dims = "group")
    k = kappa[group_idx]
    l = llambda[group_idx]
    
    # Cumulative Weibull distribution fit
    b = pm.Deterministic("b", l**-k)
    mean = pm.Deterministic("mean", 1-np.exp(-b*x**k))
    std = pm.Uniform('std',0,10)
    obs = pm.Normal('obs', mu=mean, sd = std, observed = Y)


with model:
    trace = pm.sample(1000, tune=2000, chains=2,target_accept=0.99)    

Can you post the details of the error message and provide info regarding the pymc3 version you have installed? Your code seem to be running ok for me (though I had to fake your observations dataframe).

I appreciate you taking a look at this.

I get the following error message:

ValueError: operands could not be broadcast together with shapes (61,) (1175,) 

My “group” index has 61 elements (61 observation groups) and “build” has 4 types, there are 1175 data points corresponding to “obs_id”. Obviously there would be a dimension mismatch if I were trying to do any sort of dot product, but I’m perplexed why the shape is being broadcasted as (61,) rather than 1 since I am selecting a single value from the index.

I am currently running version 3.10.0.

I think there may be a conflict in how you are requesting your normals (kappa and lambda). In each case, it looks like you may be supplying 1175 means (e.g., m1[build_idx]) but also specifying that the shape be (61,) (because dims = "group"). So maybe you intended this?

kappa = pm.Normal("kappa",
                  m1[build_idx],
                  s1[build_idx],
                  dims = "obs_id")

Alternatively, I could imagine that you are trying to create a hierarchical model with one kappa per build, each parameterized by the corresponding group (or one per group and parameterizing each by the appropriate build, I’m not sure which way your hierarchy works). So something like this:

kappa = pm.Normal("kappa",
                  m1[build_to_group_map],
                  s1[build_to_group_map],
                  dims = "build_idx")

You can’t use build_idx to do this, because it maps observations to builds.

Of course, I may be way off. Labeled coordinates and dimensions are relatively new and I am definitely still getting comfortable with them myself. @OriolAbril provided a couple of nice examples based on classic models/data sets.

1 Like

It looks like an issue with how the hierarchies are implemented. Both build_idx and group_idx have dim obs_id so both will have the same length of 1175. As your code runs, I am assuming that build_idx is an array of length 1175 whose values are integers between 0-3, therefore, when doing m1[build_idx] you are defining an array of length 1175 whose values at each position are the m1 value correponding to the index of build_idx at the same position, this is not compatible with the dimensions you have defined for kappa and so the code breaks.

I think the second solution provided above by @cluhmann will fix your issue when it comes to indexing, build_to_group_map should be an array of length 61 whose values are between 0-3 and indicate the build that corresponds to each group. It’s dims should then continue being group. And then you’ll be able to use group_idx to index assuming it is an array of length 1175 whose values go from 0 to 60.

1 Like

Thank you so much @OriolAbril and @cluhmann. You corrected a big misunderstanding on my part on the behavior of the labeled coordinates. I followed the suggestion of leveraging an intermediate variable ‘build_to_group_map’, and that did the trick. Working code below.

Y = observations['z']
X = observations['x']
build_idx, build = pd.factorize(observations["building"], sort=True)
group_idx,group = pd.factorize(observations['group'],sort=True)

build_group_map = {}
for i,g in enumerate(group_idx):
  build_group_map[g] = build_idx[i]
build_to_group_map = list(build_group_map.values())

coords = {"group":group,"build":build,"obs_id":range(observations.shape[0])}
with pm.Model(coords=coords) as model:
    
    # Setup indexing variables
    build_idx = pm.Data("build_idx", build_idx, dims="obs_id")
    group_idx = pm.Data("group_idx", group_idx, dims="obs_id")
    build_to_group_map = pm.Data("build_to_group_map", build_to_group_map, dims="group")
    x = pm.Data("x", X)
    
    # Hyper-priors
    m1 = pm.Uniform('m1',0,10, dims = 'build')
    m2 = pm.Uniform('m2',0,10, dims = 'build')
    s1 = pm.Uniform('s1',0,2, dims = 'build')
    s2 = pm.Uniform('s2',0,2, dims = 'build')
        
    kappa = pm.Normal("kappa", m1[build_to_group_map], s1[build_to_group_map], dims = "group")
    llambda = pm.Normal("llambda",m2[build_to_group_map], s2[build_to_group_map], dims = "group")
    
    # Cumulative Weibull distribution fit
    k = kappa[group_idx]
    l = llambda[group_idx]
    b = pm.Deterministic("b", l**-k)
    mean = pm.Deterministic("mean", 1-np.exp(-b*x**k))
    std = pm.Uniform('std',0,10)
    obs = pm.Normal('obs', mu=mean, sd = std, observed = Y)


with model:
    trace = pm.sample(1000, tune=2000, chains=2,target_accept=0.99)    
2 Likes

I know nothing about the model, but it may be a good idea to check that all observations with group g have the same build when looping to create the dict (unless it has already been done at some point before). You could use a groupby instead of a loop to check this and get the results in 3 steps: define groupby, check all elements of building are the same within a group and then get an array with the 1st value in each group to use as build_group_map (as they are the same it does not matter which).

Glad that helped. The coords were a red herring, but I do think that the dimensionality of parameters is slightly more implicit when using the coordinates (e.g., dims='my_dimension') than when specifying the shape explicitly (e.g., shape=len(df['group'].unique())). But that very well may be due to my relative inexperience with coords/dims. Maybe I’ll feel differently as I continue to use the new features.