[R] non-ideal behavior in princomp/ not a feature but a bug

[Prof Brian D Ripley]

On Fri, 29 Sep 2000, Ritter, Christian C SRTCL-CTGAS wrote:

> ... I checked and Brian and I are both right (see bottom for prior mail
> exchange). 

Um: you only appeared to be right because you use an obselete version of R
(and I don't think you said so), whereas I used the current release. Moral:
always upgrade before posting bugs.


> Let me explain:
> =============================================================
> 1. Indeed, in principle, princomp allows data matrices with are wider than
> high.
> Example:
> > x1
>      [,1] [,2] [,3] [,4]
> [1,]    1    1    2    2
> [2,]    1    1    2    2
> > princomp(x1)
> Call:
> princomp(x = x1)
> 
> Standard deviations:
> Comp.1 Comp.2 Comp.3 Comp.4 
>      0      0      0      0 
> 
> 4  variables and  2 observations.
> 
> So far so good
> ===========================================================
> 
> 2. However, not all matrices work, such as the first one I tried yesterday.
> Here is
> an example of a matrix which does not work.
> > x
>      [,1] [,2] [,3] [,4]
> [1,] -0.4 -1.6 -1.0 -1.8
> [2,] -2.2 -0.3  1.4  0.2
> > princomp(x)
> Error in princomp(x) : covariance matrix is not non-negative definite
> > 
> 
> Seeing this behavior on the matrix I tried, I hastily concluded that
> princomp was not
> suitable for this type (as similar functions are in many other software
> packages which
> are written for psychologists but not for chemometricians).
> 
> =======================================================
> However, the problem is not with the design of princomp in general, but with
> a little
> check it performs inside, which can have problems with roundoff:
> 
>     if (any(edc$values < 0)) 
>         stop("covariance matrix is not non-negative definite")

Use the current version of R: the check is actually

        if (any(ev[neg] < -9 * .Machine$double.eps * ev[1])) 

Now, that may not be adequate (I don't know where the precise value
comes from) and had I written it I would have used more like 100 than 9.
But it covered all the examples I tried.

> As it happens, the covariance of x1 has eigenvalues of:
> > eigen(cov.wt(x1)$cov)$values
> [1] 0 0 0 0
> and x has eigenvalues
> > eigen(cov.wt(x)$cov)$values
> [1]  7.345000e+00  1.110223e-16  0.000000e+00 -1.110223e-16
> Therefore, the observed behavior is logical (but wrong).
> 
> Oh, by the way:
>          _         
> platform Windows   
> arch     x86       
> os       Win32     
> system   x86, Win32
> status             
> major    1         
> minor    1.0       

That's not current.

[...]

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272860 (secr)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._