I have a multivariate distribution parametrized by a vector
mu and a Cholesky decomposed covariance matrix
sigma = L.L^T. Pretty much the same as in the
Now I want to update the prior whenever I have new data available, but the examples do not cover how to do this for cases where the parameters are vector valued.
I looked at
this question which discusses the problem but no clear solution is provided.
What is the right way to do this with PyMC?
Is it possible to apply
from_posterior to each axis of the trace samples and then join the Interpolated distributions? How would this join/concatenation work?
You mean the updating prior using interpolation? Yeah that wouldnt work in your case with a vector prior.
But maybe you can parameterized on the things like eta of LKJCorr and update that instead?
Oh maybe try something like:
I had a look at your code. Unfortunately, I dont see much more room to improve the speed, as the Interpolation rv always need to reinitialized, and the bottom neck in the evaluation of scipy Spline smooth function is also difficult to improve without rewriting everything in theano.
Alternatively, if you are happy with using some other distribution as an approximation, it would simplify the code a lot:
new_data = pd.DataFrame(data=[generate_data() for _ in range(1)], columns=generate_data().key…
Thank you, I think I mostly get it.
I can get the means of each set of trace samples for Xi in
X = [X1, …Xn] and then use a student-t distribution to approximate the actual distribution over X.
Assuming that my understand is correct, there are two questions:
How would I provide the n mean values to a single
How would the covariance (sd?) be modeled?
Maybe it would be easier if you can share your code and indicate which are the parameters you would like to do posterior update.
A notebook can be found
I doubt that the way I’ve done this in the notebook is correct because I wasn’t sure what the right way to approximate
Hmmm i see. In theory you can always approximate the latent freeRVs (not necessary the RV you wrote down, as pymc3 sample in the nonbounded real space) with a distribution, but it would involves changing the model quite a bit…
Thanks for your reply. Could you please elaborate a little more? I understand your comment about approximating latent free RVs but do you have any ideas on how I could go about modeling that?
@junpenglao I re-parametrized my model, like so:
with pm.Model() as model1:
mu = pm.InverseGamma('mu', 3., 1., shape=shape)
eta = pm.Gamma('eta', 3., 1.)
sd_dist = pm.HalfCauchy('sd_dist', 2.5)
packed_L = pm.LKJCholeskyCov('packed_L', n=shape, eta=eta, sd_dist=sd_dist)
L = pm.expand_packed_triangular(shape, packed_L)
sigma = pm.Deterministic('sigma', L.dot(L.T))
pm.MvNormal('observation', mu, sigma, observed=observed)
But this gives the following error:
~/.local/lib/python3.5/site-packages/pymc3/distributions/multivariate.py in logp(self, x)
981 sd_vals = tt.sqrt(variance)
--> 983 logp_sd = self.sd_dist.logp(sd_vals).sum()
984 corr_diag = x[diag_idxs] / sd_vals
AttributeError: 'TransformedRV' object has no attribute 'logp'
It seems that
sd_dist must always be a
pm.Distribution, wouldn’t it be more useful if I could model it like all the other parameters?
Sorry about the non response. we did not directly model the var parameters on the diag - but maybe you can model the diag and the corr matrix separately.
Thanks, is there a reason for this choice of design?