[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Multivariate pdf of a normal distribution
From: |
Mike Miller |
Subject: |
Re: Multivariate pdf of a normal distribution |
Date: |
Sun, 6 Nov 2005 19:16:41 -0600 (CST) |
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
-------------------------------------------------------------
- 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