Draw number of MvNormal in once

Questions
1

Dear pymc3 community.

Hello, I’m currently working discrete choice modelling with pymc3.

I built several basic models to build up my model.

So, the point is that I need to draw the number of MvNormal coefficients in pooled hierarchical model

I already build hierarchical bayesian logit model with following methods.

with pm.Model() as modellogit:
    # beta coefficient variance  
    packed_L = pm.LKJCholeskyCov('packed_L', 
                                 n= beta_size, 
                                 eta = 2, 
                                 sd_dist = pm.HalfCauchy.dist(1.5))
    Lower = pm.expand_packed_triangular(beta_size,packed_L)
    Sigma = pm.Deterministic('Sigma',pm.math.dot(Lower,Lower.T))
    
    # Beta coefficient
    betabar = pm.MvNormal('betabar',mu = np.zeros(beta_size),
                          cov =np.identity(beta_size),
                          shape=(1,beta_size))
    
    beta = pm.MvNormal('beta',mu = betabar,
                       cov = Sigma,
                       shape = (nid,beta_size))
    
    temp1 = tt.sum(attr*beta[df1['id'].values,],axis=1)
    temp2 = tt.reshape(temp1,(nid*ntrial,nalt))

However, the thing I confronted is that each individual shold get different covariance value.

So I replace
Sigma to (nid,beta_size,beta_size) tensor3 objects and
betabar to (nid,beta_size) matrix objects

to assign different sigmas per individual.

and tried it. but it does not work due to cov needs to 2 dim matrix

And then. I tried for loop in below way.

beta = []
for i in range(nid):
    Sigma = (... new Sigma...)
    _beta = pm.MvNormal('beta %d' %i ,mu = betabar,
                       cov = Sigma,
                       shape = (1,beta_size))
    beta.append(_beta)
beta = tt.stack(beta)
print(beta)
    
temp1 = tt.sum(attr*beta[df1['id'].values,:],axis=1)
temp2 = tt.reshape(temp1,(nid*ntrial,nalt))

But it also did not work since temp1 part produces (nid * ntrial * nalt * beta_size) array.
(attr matrix is size of (nid * ntrial * nalt, beta_size ) array)

How could I solve this problem?

I’m very appreciate your contribution to bayesian world!

2

I see two MvNormal

# Beta coefficient
    betabar = pm.MvNormal('betabar',mu = np.zeros(beta_size),
                          cov =np.identity(beta_size),
                          shape=(1,beta_size))
    
    beta = pm.MvNormal('beta',mu = betabar,
                       cov = Sigma,
                       shape = (nid,beta_size))

Which one has the problem?

3

the code first written is worked well. even though it is slow.

So I changed beta part of first code into below loop way.

If I stack beta into one matrix, then the temp1 part of below code does not work at all.