I’m trying to fit a multivariate normal analogous to this example
Different from the example, I have several sets of observations for which I want to estimate independent parameters. I can’t seem to figure out how to add a dimension to the example to handle this (and I want to avoid an ugly loop thing). Normally I know how to work with the shapes but the observations being pairs somehow makes it different…
I’ve tried this (among other things)
t, trials = pd.factorize(baseline_data['Trial Index'], sort=True) with pm.Model(coords=coords) as model: chol, corr, stds = pm.LKJCholeskyCov('chol', n=2, eta=2., sd_dist=pm.HalfCauchy.dist(500,shape=(len(trials), 2)), compute_corr=True) μ = pm.Normal('μ', 0., 2000., testval=baseline_data['X'].mean(axis=0), shape=(len(trials), 2)) obs = pm.MvNormal('obs', μ[t], chol=chol[t], observed=np.array([baseline_data_static['X'].values, baseline_data_static['Y'].values]).T) trace = pm.sample(tune=1000, return_inferencedata=True, init='advi+adapt_diag')
It results in "
IndexError: index 2 is out of bounds for axis 0 with size 2"
I have also tried different orders of the dimensions. And somewhere there is a version of the example with named dimensions but that didn’t help either…
Any suggestions are highly appreciated!