4:28 p.m.
Hi,
perhaps there's a bug in princomp in the case of reduced rank. See this
illustration:
<hansl>
matrix A = { 0, 0, 0, 0; \
0.0000 , 2.2790, -1.5619, -1.6955; \
0.0000 , -1.5619 , 1.2914 , 1.6958 ; \
0.0000 , -1.6955 , 1.6958 , 2.6103 }
eval mcov(A) # OK
eval mcorr(A) # OK (correct NAs)
matrix v = {}
eval eigensym(A, &v) # OK
print v
eval princomp(A, 1, 1) # try cov; only gives NAs
eval princomp(A, 1) # fails
</hansl>
In my understanding a PC is basically just an eigenvector of the
covariance (or correlation) matrix, so if eigensym(A) is able to
calculate some eigenvectors for the non-zero eigenvalues, so should be
princomp(), no?
BTW, the error message for the last line is "syntax error", which is
also wrong I think.
Thanks,
sven
6:10 p.m.
On Wed, 26 Jul 2017, Sven Schreiber wrote:
Hi,
perhaps there's a bug in princomp in the case of reduced rank. See this
illustration:
<hansl>
matrix A = { 0, 0, 0, 0; \
0.0000 , 2.2790, -1.5619, -1.6955; \
0.0000 , -1.5619 , 1.2914 , 1.6958 ; \
0.0000 , -1.6955 , 1.6958 , 2.6103 }
eval mcov(A) # OK
eval mcorr(A) # OK (correct NAs)
matrix v = {}
eval eigensym(A, &v) # OK
print v
eval princomp(A, 1, 1) # try cov; only gives NAs
eval princomp(A, 1) # fails
</hansl>
In my understanding a PC is basically just an eigenvector of the covariance
(or correlation) matrix, so if eigensym(A) is able to calculate some
eigenvectors for the non-zero eigenvalues, so should be princomp(), no?
Not sure, I'll defer to Jack on that. But I note that R's princomp()
and prcomp() both choke on these data ("must be finite").
BTW, the error message for the last line is "syntax error",
which is also
wrong I think.
Looks like it, yes. That should be fixable.
Allin
10:41 a.m.
On Wed, 26 Jul 2017, Allin Cottrell wrote:
On Wed, 26 Jul 2017, Sven Schreiber wrote:
> Hi,
> perhaps there's a bug in princomp in the case of reduced rank. See this
> illustration:
[...]
Not sure, I'll defer to Jack on that. But I note that R's
princomp() and
prcomp() both choke on these data ("must be finite").
I guess that your script is not using princomp() in the intended way: what
you're supplying the function is the _covariance_ matrix of the data, not
the _data_; if you actually had a dataset with 4 variables and 4
observations, and variable 1 is all zeros, then it would have zero
variance, so I guess that in any real-life cases it wouldn't be considered
in the first place. But if you insist in doing so, the first principal
component on the covariances can be returned ok, as
<hansl>
matrix A = { 0, 0, 0, 0; \
0.0000 , 2.2790, -1.5619, -1.6955; \
0.0000 , -1.5619 , 1.2914 , 1.6958 ; \
0.0000 , -1.6955 , 1.6958 , 2.6103 }
matrix v = {}
matrix S = mcov(A)
eval eigensym(S, &v)
print v
eval A * v[,4]
</hansl>
OTOH, in this case a correlation matrix cannot be meaningfully defined, so
it's only fair that princomp() should choke.
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
11 a.m.
Two more things on this; one, I wrote
But if you insist in doing so, the first principal component on the
covariances can be returned ok, as
<hansl>
matrix A = { 0, 0, 0, 0; \
0.0000 , 2.2790, -1.5619, -1.6955; \
0.0000 , -1.5619 , 1.2914 , 1.6958 ; \
0.0000 , -1.6955 , 1.6958 , 2.6103 }
matrix v = {}
matrix S = mcov(A)
eval eigensym(S, &v)
print v
eval A * v[,4]
</hansl>
which is wrong, because you should centre the data first, and use
cdemean(A) in place of A.
Two, I found that princomp actually uses standardised data when computing
its output even when the third argument is non-zero, which is IMO wrong.
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------
11:17 a.m.
Am 28.07.2017 um 13:00 schrieb Riccardo (Jack) Lucchetti:
Two more things on this; one, I wrote
> But if you insist in doing so,
I know that princomp() should be applied to data, not to the cov-matrix.
This is just a minimal numerical example to showcase the issue (and with
some path dependence also playing a role).
the first principal component on the
> covariances can be returned ok,
Yes that's exactly what I meant and why I didn't understand that
princomp() doesn't manage to do that.
which is wrong, because you should centre the data first, and use
cdemean(A) in place of A.
OK, but AFAICS doesn't change anything of the princomp() problems here.
Two, I found that princomp actually uses standardised data when
computing its output even when the third argument is non-zero, which is
IMO wrong.
Could be. I vaguely remember some of those things were discussed when
the third optional argument was introduced (the covariance-based variant).
cheers,
sven
2708
days inactive
2710
days old
4 comments
3 participants
participants (3)
-
Allin Cottrell -
Riccardo (Jack) Lucchetti -
Sven Schreiber