Thank you! This works great for one of my use cases.
However, I used it in a bit of an unconventional way. I noticed that the way you used the MvNormal is not too different from the prior_from_idata implementation, so it may be possible to use it this way:
mu_m, mu_s = sp.stats.norm.fit(az.extract(idatas[i].posterior)['mu'].values)
sigma_m, sigma_s = sp.stats.norm.fit(az.extract(idatas[i].posterior)['sigma'].values)
prior = prior_from_idata(idatas[i], var_names=['mu', 'sigma'])
mu = prior['mu'] * mu_s + mu_m
sigma = prior['sigma'] * sigma_s + sigma_m
As the prior_from_idata would apply the Mvnormal approx. to the whole joint posterior and then we can extract it per each singular posterior (i.e. prior[‘mu’], prior[‘sigma’]). Then, the mean and SD obtained from sp.stats.norm.fit, as in your approach, could be directly applied to the prior_from_idata output. Please correct me if I got this wrong. Maybe I got it working just by chance.