I have a matrix representing the distance between random normal variables. I’d like to estimate an additional variable representing the rate of decline in co-variance between these variables based on their distance. The centered implementation looks like this:
distance_matrix = np.array([[0,1,2],
[1,0,1],
[2,1,0]])
# centered implementation
rho_sq = pm.HalfNormal('rho_sq', sd=1) # rate of decline by distance
cov = pm.math.exp(-rho_sq * distance_matrix**2)
pm.MvNormal('correlated_vars', mu=0, cov=cov, shape=len(distance_matrix))
The non centered implementation requires a cholesky decomposition. How can I do that in pymc3?
distance_matrix = np.array([[0,1,2],
[1,0,1],
[2,1,0]])
# non centered implementation
rho_sq = pm.HalfNormal('rho_sq', sd=1) # rate of decline by distance
cov = pm.math.exp(-rho_sq * distance_matrix**2)
# this does not work
chol = np.linalg.cholesky(cov)
uncorrelated_vars = pm.Normal('uncorrelated_vars', mu=0, sd=1, shape=len(distance_matrix))
correlated_vars = pm.Deterministic('correlated_vars', tt.dot(chol, uncorrelated_vars.T))