Any suggestions on continous error "ValueError: shape mismatch: value array of shape (62,) could not be broadcast to indexing result of shape (0,)"

9

I put the revised code for my notes

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pymc as pm
import pytensor.tensor as pt
import arviz as az
import pytensor
import time

D_phylo = pd.read_csv('Combo_Dmatrix.csv')
COMBO_OTU = pd.read_csv('COMBO_OTU.csv')
y_obs = pd.read_csv('ynew3.csv')
y_true = pd.read_csv('truevalues.csv')

D_phylo = D_phylo.to_numpy()
D_phylo = D_phylo[:,1:]
COMBO_OTU = COMBO_OTU.to_numpy()
COMBO_OTU = COMBO_OTU[:,1:]

y_obs = y_obs.to_numpy()
y_obs = y_obs[:,1:]
y_true = y_true.to_numpy()

K1_combo = int(62 * 0.85)
# rho = 0.25
# sigma = 1
cov =  1**2* np.exp(-2* 0.25* D_phylo)
samples = y_obs.shape[0]
specise = y_obs.shape[1]
phai = np.zeros([samples,1])

N = np.sum(COMBO_OTU[0:samples,:], axis=1).reshape(-1, 1)

tt = samples
chain_length = 2000
tune_length = 1000
spe = ['spe{}'.format(i) for i in range(specise)]
sam = ['sam{}'.format(i) for i in range(tt)]

coords = {"specie": spe, "sample": sam}
ind = np.arange(specise)
ind = pt.as_tensor_variable(ind)
cov = pt.as_tensor_variable(cov)

pp = pt.as_tensor_variable(np.zeros([tt,specise]))
true_abundance = pt.as_tensor_variable(np.zeros([tt,specise]))
obs_like = np.array(y_obs[:,:].tolist())
NN = np.sum(np.array(y_obs[:,:].tolist()),axis = 1)



zk = np.arange(specise, dtype='int32')
zk = pt.as_tensor_variable(zk)
zp = np.ones([specise], dtype='int64')*int(62)
zp = pt.as_tensor_variable(zp)

start = time.time() 


def def_test2(i,pll, U,expend_con,ind_tem):
    p2 = pll
    tem = pm.math.sum(pm.math.exp(U[i, :p2-1]))
    concatenated_values = pm.math.concatenate((pm.math.exp(U[i, :p2-1]) / (1 + tem), [1 / (1 + tem)]))
    expend_con = pytensor.ifelse.ifelse(pt.eq(p2,0),
                                        pt.set_subtensor(expend_con[i,:], pt.zeros(62)+1/62),
                                        pt.set_subtensor(expend_con[i, ind_tem], concatenated_values))
    return expend_con
  
if __name__ == '__main__':
    with pm.Model(coords=coords) as Multi:        
        #########Prior########
        #alpha = pm.Normal('alpha',mu =  5,   sigma = 1,  shape = specise)       
        theta = pm.Normal('theta', mu = 1 , sigma = 1, shape = specise)
        
        # T = pm.Gamma('T',7.5,1) #15  gamma or lognormal
        T = pm.Gamma('T',0.01,0.01) #15  gamma or lognormal
        #D = pm.HalfNormal('D',sigma = 1, shape = specise)
        p = pm.Deterministic('p',1-pm.math.exp(-D_phylo[-1,:]/T)) # the last column is the reference taxon
        delta = pm.Bernoulli("delta",p = p, shape = specise)   #np.ones([specise])
        # print(time.time()-start)          
        
        ind_tem = ind[pm.math.lt(delta[:],1)]
        pll = ind_tem.shape[0] 
        
        z_ind_tem = zk   
        ind_tem = pytensor.ifelse.ifelse(~pt.all(pll), z_ind_tem, ind_tem)
        zpll =    pytensor.ifelse.ifelse(~pt.all(pll), ind_tem.shape[0], pll)
        
        th_t = []
        ct = []
        for i in range(1):
            theta_tem = pt.as_tensor_variable(np.ones([specise]))
            cov_tem = pt.as_tensor_variable(np.diag(np.ones([specise])))
            t1 =pt.set_subtensor(theta_tem[0:zpll], theta[ind_tem])
            t2 = pt.set_subtensor(cov_tem[0:zpll,0:zpll], cov[ind_tem][:, ind_tem])
            th_t.append(t1) 
            ct.append(t2)  
               
        U = pm.MvNormal("U", mu = th_t, cov = ct, dims = ("sample", "specie"))            
                               
        i = pt.arange(tt)
        expend_con = pt.zeros([tt,specise])+1e-10
        Results, _ = pytensor.scan(def_test2,
                                    sequences = [i],
                                    non_sequences = [pll, U,expend_con,ind_tem])
        
        for i in range(tt):
            expend_con = pt.set_subtensor(expend_con[i,:],Results[i,i,:])
              
        ########Likelihood######## 
                         
        likelihood_y = pm.Multinomial('likelihood', n =NN[0:tt], p = expend_con, observed = obs_like[0:tt,:]) #N should be total counts of species of each sample
        print(time.time()-start) 
        ########sampling######## 
        trace = pm.sample(chains = 4, cores = 8,  step = pm.Metropolis()) #chain_length, tune = tune_length,
        az.plot_trace(trace)
        
        ########output########
        summary = az.summary(trace,hdi_prob=0.95,stat_focus = 'mean')
        print(summary)