Sturla Molden <sturla <at> molden.no> writes: > > Yes, this is what I am computing. I am computing the pdf of a very high- > > dimensional multivariate normal. Is there a specialized method to compute > > this? > > If you use cho_solve and cho_factor from scipy.linalg, you can proceed > like this: > > cx = X - m > sqmahal = (cx*cho_solve(cho_factor(S),cx.T).T).sum(axis=1) > > where X is the data (n x p), m is mean (1 x p) and S is the covariance (p > x p). > > We could do this twice as fast using a hypothetical method tri_solve that > calls LAPACK subroutine DTRTRS.
Sorry for the laziness of this question... As a reminder, I am working with high-dimensional data (10K dimensions). I was computing the log of the MVN pdf because the probabilities would almost all go to zero. Do you suppose the method you have shown me will be numerically stable (the probabilities will be small but they stay above zero)? By the way, I would be happy to implement the method you desire in a couple months if you are willing to do a little hand-holding. Thanks for the excellent suggestion. - dan _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
