Hierachi sample error

Questions version agnostic
1

hi,thank u in advance

import pymc as pm
import numpy as np
import pandas as pd
import arviz as az
import matplotlib.pyplot as plt
import xarray as xr
import aesara

#data
trainSize,testSize = 1000,400 
group_num_all = 7
group_num_train = 5
field_num = 3
group_size = (trainSize+testSize)/group_num_all
group_idx_val = np.repeat(range(group_num_all),group_size)
group_idx_train,group_idx_test = group_idx_val[:trainSize],group_idx_val[trainSize:]
x_train,x_test=np.random.random((trainSize,field_num)),np.random.random((testSize,field_num))
y_train,y_test=np.random.random((trainSize,1)),np.random.random((testSize,1))

#modeling 

with pm.Model() as model:
    slope_mu = pm.Normal("slope_mu", 0.5, 0.1)
    slope_sigma = pm.Normal("slope_sigma", 0.5, 0.1)
    intercept_mu = pm.Normal("intercept_mu", 0.5, 0.1)
    intercept_sigma = pm.Normal("intercept_sigma", 0.5, 0.1)
    
    group_num = pm.MutableData("group_num",group_num_train)
    slope = pm.Normal("slope", slope_mu, slope_sigma,shape=(group_num,field_num))
   
   #note key is here---,
   # if shape=(group_num,) then sample is ok,but sample_prior_predictive fails
   # if shape=(group_num,1) then sample fails,but sample_prior_predictive is ok
    intercept = pm.Normal("intercept", intercept_mu, intercept_sigma,shape=(group_num,1))
    
    group_idx = pm.MutableData("group_idx",group_idx_train)
    x = pm.MutableData("x",x_train)
    mu = intercept[group_idx] + pm.math.sum(x*slope[group_idx],axis=1)
    
    sigma = pm.Exponential("sigma", 1.0)
    pm.Normal("obs", mu=pm.Deterministic("y",mu), sigma=sigma, observed=y_train)

    idata = pm.sample_prior_predictive(random_seed=100)
    trace = pm.sample()

sample error:
ValueError: array is not broadcastable to correct shape
Apply node that caused the error: AdvancedIncSubtensor1{inplace,inc}(Elemwise{Composite{((-i0) / i1)}}.0, Elemwise{true_div,no_inplace}.0, group_idx)
Toposort index: 77
Inputs types: [TensorType(float64, (None, None)), TensorType(float64, (1000, None)), TensorType(int32, (None,))]
Inputs shapes: [(5, 1), (1000, 1000), (1000,)]
Inputs strides: [(8, 8), (8000, 8), (4,)]

sample_prior_predictive error:
ValueError: Output size (1000, 1) is not compatible with broadcast dimensions of inputs (1000, 1000).

2

Your example ran fine for me on Colab after I reduced the train/test sizes to 10 / 4 from 1000 / 400 for debugging purposes. You can see that notebook here.

I think there might be an issue in the setup with group_idx_train,group_idx_test and the values in it generated from np.repeat. Could you talk a little more about the desired model structure?

3

thanks,ck. actually,I just want to test a hierrachi case with multiple dim/features,and no finding a workable example.so I tried to make one.And as you can see the log,it reports some issues about shape. I think the root cause is there is some wrong setting for shape somewhere.

4

I took a second pass at your code and attempted to “fix” it - take a look and see if it matches what you were trying to get: 2nd attempt at fixing hierarchical model example · GitHub. Thanks for providing a running end to end example. It makes it much easier to debug :slight_smile:

I have the following comments:

  • You were using group_num as a model variable to hold the number of active groups in the training data, it appears. I understand 100% why you did this but it’s unnecessary - what I would do instead is to instantiate all the groups’ intercepts/slopes in the model and then run the sampler with the training data you do have. It’s okay if there are groups without any training data in the model - PyMC and the built in samplers can handle this just fine!
  • I changed your index set variable group_idx_val from something like [0, 0, 0,...., 1, 1, 1,...., 6, 6] to a random sample of length totalSize from the set {0,...,6} to make it more reflective of how these types of models often have their data organized; if that’s not the case.
5

thanks ck.Eventually, I fixed the shape issue by “flatten”.
pm.Normal(“obs”, mu=pm.Deterministic(“y”,mu), sigma=sigma, observed=y_train.flatten()).
I have to say using pymc is really “art” :slight_smile: