Thanks for the reply. David Warde-Farley <dwf <at> cs.toronto.edu> writes: > On 2010-08-30, at 11:28 AM, Daniel Elliott wrote: > > Large matrices (e.g. 10K x 10K) > > > Is there a function for performing the inverse or even the pdf of a > > multinomial normal in these situations as well? > > There's a function for the inverse, but you almost never want to use it, especially if your goal is the > multivariate normal density. A basic explanation of why is available here: > http://www.johndcook.com/blog/2010/01/19/dont-invert-that-matrix/ > > In the case of the multivariate normal density the covariance is assumed to > be positive definite, and thus a > Cholesky decomposition is appropriate. scipy.linalg.solve() (NOT numpy.linalg.solve()) with the > sym_pos=True argument will do this for you.
You don't think this will choke on a large (e.g. 10K x 10K) covariance matrix? Given what you know about how it computes the log determinant and how the Cholesky decomposition, do you suppose I would be better off using eigen- decomposition to do this since I will also get the determinant from the sum of the logs of the eigenvalues? > What do you mean by a "multinomial normal"? Do you mean multivariate normal? Offhand I can't remember if it > exists in scipy.stats, but I know there's code for it in PyMC. Oops, I meant multivariate normal. I found code in the pymix library and so far I like what I see there. I am sure it works well for most cases. However, with the exception of one R library, I have never found a packaged implementation that works on covariance matrices as large as what I am talking about above. _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
