I am trying to construct an inference model where the N observables of multivariate normal with M dimension. Each data point comes with its own Error covariance, cov_err which we know, and an intrinsic covariance, cov_int which we want to learn.
cov = np.array([cov_int + cov_err[i] for i in range(N)]) yObs = pm.MvNormal('yObs', mu=mu, cov=cov, observed=obs)
but MvNormal does not support NxMxM type covariance. I was wondering if there is a fast solution for this problem assuming N ~ 100k, M ~ 10.