Here is the (not working) model to give you a better idea of what I’m trying to do. ‘invP’ is the inverse of a known covariance matrix ‘sigmaP’. ‘invL’ is the inverse of the covariance matrix (one for each participant) that I’m trying to estimate. ‘muij’ is the center of a distribution with variability ‘tau^2’, that is computed using these inverse matrices (invP and invL). ‘targets’ is a n_subjects x n_trials(=60) by 2 array containing (observed) x,y values that were sampled from a multivariate normal with mean=‘muP’ and cov=‘sigmaP’. The subject has to memorize the target’s location and reproduce it from memory. ‘endpts’ is a n_subjects x n_trials by 2 array containing the observed responses that I’m fitting to.
The first issue I don’t know how to tackle is to compute ‘muij’ for each participant and trial. I’m getting an error at this step. Also, as per my previous reply, I still don’t understand how this model is computing a separate ‘chol’ (and ultimately ‘invL’) for each participant. It looks to me like only one covariance matrix is estimated.
with pm.Model() as PiN:
sd_dist = pm.InverseGamma.dist(alpha=1., beta=1)
chol, corr_, stds_ = pm.LKJCholeskyCov('lkjcc', n=2, eta=2., sd_dist=sd_dist, compute_corr=True)
cov = pm.Deterministic('cov', chol.dot(chol.T))
invP = pm.Deterministic('invP', tt.as_tensor_variable(np.linalg.inv(sigmaP)))
invL = pm.Deterministic('invL', pm.math.matrix_inverse(chol))
μij = pm.Deterministic('μij', pm.math.dot(pm.math.matrix_inverse(invP + invL), pm.math.dot(invP, muP) + pm.math.dot(invL, tt.as_tensor_variable(targets[subs_idx]))))
τ = pm.HalfNormal('τ', sd=1)
τ_cov = pm.Deterministic('τ_cov', tt.as_tensor_variable(τ**2 * tt.as_tensor_variable(np.identity(2))))
mvn_dist = pm.MvNormal.Dist(mu=uij, chol=τ_cov, shape=(nsubs, 2))
sij = pm.Potential('sij', (mvn_dis.logp(endpts[1,subs_idx])))