[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Multivariate pdf of a normal distribution
From: |
Henry F. Mollet |
Subject: |
Re: Multivariate pdf of a normal distribution |
Date: |
Sun, 06 Nov 2005 19:35:33 -0800 |
User-agent: |
Microsoft-Entourage/11.1.0.040913 |
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
-------------------------------------------------------------
- Re: Multivariate pdf of a normal distribution, (continued)
Re: Multivariate pdf of a normal distribution, Michael Creel, 2005/11/07
Re: Multivariate pdf of a normal distribution, Gorazd Brumen, 2005/11/06