I am trying to duplicate the hierarchical model shown above in Pymc3.
I have no problem to build this model in pymc3 with the residuals for each OS and PD follow by N~(0,sigma^2). After reading this Analysis of AR(1) in Pymc3, I still do not know how to incorporate the first order Autoregressive residuals AR(1) into my model in Pymc3.
Therefore, I have no choice but to post this question here.
Below is the model I built without the AR(1) as the residuals and feel free to give me advice or feedback on my model as well. Thank you for your precious time! Appreciate it!
mu_vector_prior = np.array([0.9,0.5]) cov_matrix_prior = np.array([[0.1,0],[0,0.1]]) mu_vector = np.array([5,0.5]) cov_matrix= np.array([[0.5,0],[0,0.1]]) with pm.Model() as hierarchical_model: #Prior for HyperParameters prior_mu = pm.MvNormal('prior', mu = mu_vector_prior , cov = cov_matrix_prior, shape = 2) prior_cov = pm.LKJCholeskyCov('packed_L', n=2, eta=2., sd_dist=pm.HalfCauchy.dist(2.5)) L = pm.expand_packed_triangular(2, prior_cov) sig = pm.Deterministic('sig', L.dot(L.T)) #Model's parameter #random-effect model_param_rnd = pm.MvNormal("rnd_var", mu = prior_mu, chol = L, shape = (10,2)) #fixed-effect model_param_fix = pm.MvNormal("fix_var", mu = mu_vector, cov = cov_matrix, shape = 2) #Model error epsOS = pm.Uniform('epsOS', upper = 10000, lower = 0) epsPD = pm.Uniform('epsPD', upper = 5000, lower = 0) #Expected value RLR = model_param_rnd[npcohort,0] RRF = model_param_rnd[npcohort,1] Ker = model_param_fix Kp = model_param_fix A = (P*Ker*RLR)/(Ker-Kp) B = (P*Ker*RLR*RRF)/(Ker-Kp) OS = A*((np.exp(-Kp*nptime))-np.exp(-(Ker)*nptime)) PD = B*((Ker*(1-np.exp(-Kp*nptime))) - (Kp*(1-np.exp(-Ker*nptime)))) #likelihood y_OS = pm.Normal('y_OS', mu = OS, sd = epsOS, observed = dataOS.Loss.values) y_PD = pm.Normal('y_PD', mu = PD, sd = epsPD, observed = dataPD.Loss.values)