Re: Multivariate pdf of a normal distribution

[Henry F. Mollet]

octave:27> [E,L] = eig(ones(4))
E =

   -1.7253e-01    4.6929e-01    7.0711e-01    5.0000e-01
   -1.7253e-01    4.6929e-01   -7.0711e-01    5.0000e-01
   -4.9115e-01   -7.1328e-01   -6.8934e-17    5.0000e-01
    8.3620e-01   -2.2530e-01    9.7203e-18    5.0000e-01

L =

  -0.00000   0.00000   0.00000   0.00000
   0.00000  -0.00000   0.00000   0.00000
   0.00000   0.00000  -0.00000   0.00000
   0.00000   0.00000   0.00000   4.00000

GNU Octave, version 2.1.71 (powerpc-apple-darwin8.1.0).
Why do I not get the same results?
octave:30> det(A)
ans = 0
octave:31> rank(A)
ans = 1
Only look at non-zero eigenvalue and corresponding eigenvector?
Henry


on 11/6/05 5:16 PM, Mike Miller at address@hidden wrote:

> On Sun, 6 Nov 2005, Andy Adler wrote:
> 
>> On Sun, 6 Nov 2005, Mike Miller wrote:
>> 
>>> When inv(A) exists,
>>> 
>>> inv(A') = inv(A)'
>>> 
>>> because the inverse equals the transpose of the adjoint divided by the
>>> determinant.  The adjoint of the transpose equals the transpose of the
>>> adjoint and the determinant is the same for A and A'.  The identity really
>>> follows from the fact that the determinant is not affected by the
>>> transpose operation -- the adjoint is a collection of determinants.
>>> 
>>> I agree that this must be in most matrix algebra books.
>> 
>> Yes, that is true in theory. But the error is to assume that all
>> mathematical identities remain intact when done on a finite
>> precision computer.
>> 
>> Try it:
>> A= rand(2); inv(A')- inv(A)'
>> ans =
>> 
>>    7.1054e-15   -7.1054e-15
>>    8.8818e-16    0.0000e+00
> 
> 
> That wasn't an error that anyone here was making.  The original question
> showed an example where there was a small difference.
> 
> Related to this, we also were discussing recently the problem with eigen
> decomposition of singular matrices.  For example:
> 
> 
> octave:6> [E,L] = eig(ones(4))
> E =
> 
>     -3.2026e-01    7.0711e-01    3.8397e-01    5.0000e-01
>     -3.2026e-01   -7.0711e-01    3.8397e-01    5.0000e-01
>     -2.2276e-01   -5.9234e-17   -8.3689e-01    5.0000e-01
>      8.6328e-01   -4.4132e-19    6.8941e-02    5.0000e-01
> 
> L =
> 
>    -0.00000   0.00000   0.00000   0.00000
>     0.00000   0.00000   0.00000   0.00000
>     0.00000   0.00000   0.00000   0.00000
>     0.00000   0.00000   0.00000   4.00000
> 
> octave:7> A = E*sqrt(L)
> A =
> 
>     0.00000 - 0.00000i   0.00000 + 0.00000i   0.00000 + 0.00000i   1.00000 +
> 0.00000i
>     0.00000 - 0.00000i  -0.00000 + 0.00000i   0.00000 + 0.00000i   1.00000 +
> 0.00000i
>     0.00000 - 0.00000i  -0.00000 + 0.00000i  -0.00000 + 0.00000i   1.00000 +
> 0.00000i
>     0.00000 + 0.00000i  -0.00000 + 0.00000i   0.00000 + 0.00000i   1.00000 +
> 0.00000i
> 
> octave:8> A*A'
> ans =
> 
>    1.0000  1.0000  1.0000  1.0000
>    1.0000  1.0000  1.0000  1.0000
>    1.0000  1.0000  1.0000  1.0000
>    1.0000  1.0000  1.0000  1.0000
> 
> The problem here is that numerical imprecision causes A to be complex when
> we want it to be real.  If the matrix really is singular (as this in the
> example above), we can use real() to get a good answer:
> 
> octave:9> B=real(A)
> B =
> 
>     0.00000   0.00000   0.00000   1.00000
>     0.00000  -0.00000   0.00000   1.00000
>     0.00000  -0.00000  -0.00000   1.00000
>     0.00000  -0.00000   0.00000   1.00000
> 
> octave:10> B*B'
> ans =
> 
>    1  1  1  1
>    1  1  1  1
>    1  1  1  1
>    1  1  1  1
> 
> But for a function that always works, we'll have to set some tolerance
> limit and have the function fail when the smallest eigenvalue [min(L)
> above] is less than some value.  What is a good choice of a value?  Maybe
> -1.0e-13?  Anyone?  How far off can we be when the precise min(L) is zero?
> 
> Mike
> 
> 
> 
> -------------------------------------------------------------
> Octave is freely available under the terms of the GNU GPL.
> 
> Octave's home on the web:  http://www.octave.org
> How to fund new projects:  http://www.octave.org/funding.html
> Subscription information:  http://www.octave.org/archive.html
> -------------------------------------------------------------
> 




-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------