[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