[SciPy-User] struggling and fighting with linalg
josef.pktd@gmai...
josef.pktd@gmai...
Thu Oct 21 15:16:42 CDT 2010
On Thu, Oct 21, 2010 at 1:02 PM, <josef.pktd@gmail.com> wrote:
> On Thu, Oct 21, 2010 at 12:46 PM, Warren Weckesser
> <warren.weckesser@enthought.com> wrote:
>> Hi Josef,
>>
>> On Thu, Oct 21, 2010 at 11:03 AM, <josef.pktd@gmail.com> wrote:
>>>
>>> There are a lot of recommendations in the mailing list and recipes in
>>> various packages how to use linalg more efficiently.
>>> However, yesterday I wasted several hours how to do the multivariate
>>> normal distribution with cholesky factor and solve.
>>>
>>
>> What do you mean by "*do* the multivariate normal distribution" (emphasis
>> mine)? numpy.random has multivariate_normal for generating samples; I wrote
>> a little cookbook example
>> (http://www.scipy.org/Cookbook/CorrelatedRandomSamples) before I realized
>> multivariate_normal existed.
>>
>> But I think you want more than just generating samples, right?
>
> yes, main things I need for statsmodels are loglikelihood and linear
> transformation, (and separately in other application non-linear
> transformation), plus linear transformation that comes from the
> decomposition of the inverse covariance matrix (for whitening or
> standardization). The last part is where I got stuck.
>
> The extras will be conditional, marginal distributions, cdf,
> confidence regions, ... Part of these just needs to be collected (and
> tested).
more cholesky woes: upper, lower, transpose or flip ? I still haven't
figured it out and the test don't work yet.
very simple example: autoregressive process order 1, timeseries has
increasing index with time. But cholesky seems to run backwards in
time.
>>>
>>> nobs = 5
>>> autocov = 0.8**np.arange(nobs)
>>> sigma = linalg.toeplitz(autocov)
>>> sigmainv = linalg.inv(sigma)
>>>
>>> c = linalg.cholesky(sigma, lower=True)
>>> ci = linalg.cholesky(sigmainv, lower=True)
>>>
>>> print sigma
[[ 1. 0.8 0.64 0.512 0.4096]
[ 0.8 1. 0.8 0.64 0.512 ]
[ 0.64 0.8 1. 0.8 0.64 ]
[ 0.512 0.64 0.8 1. 0.8 ]
[ 0.4096 0.512 0.64 0.8 1. ]]
>>> print tiny2zero(ci/ci.max())
[[ 1. 0. 0. 0. 0. ]
[-0.8 1. 0. 0. 0. ]
[ 0. -0.8 1. 0. 0. ]
[ 0. 0. -0.8 1. 0. ]
[ 0. 0. 0. -0.8 0.6]]
>>>
>>> "this is the text book transformation for efficient GLS and MLE"
>>> print 'coefficient for first observation', np.sqrt(1-autocov[1]**2)
coefficient for first observation 0.6
>>> ci2 = ci[::-1,::-1].T
>>> print tiny2zero(ci2/ci2.max())
[[ 0.6 0. 0. 0. 0. ]
[-0.8 1. 0. 0. 0. ]
[ 0. -0.8 1. 0. 0. ]
[ 0. 0. -0.8 1. 0. ]
[ 0. 0. 0. -0.8 1. ]]
>>>
>>> print np.dot(ci/ci.max(), np.ones(nobs)) #wrong
[ 1. 0.2 0.2 0.2 -0.2]
>>>
>>> print np.dot(ci2/ci2.max(), np.ones(nobs)) #correct
[ 0.6 0.2 0.2 0.2 0.2]
Are there any rules to understand which direction cholesky goes?
Josef
>
> Josef
>
>>
>> Warren
>>
>>
>>
>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>>
>
More information about the SciPy-User
mailing list