SCM

SCM Repository

[matrix] View of /pkg/tests/poMatrix.Rout.save
ViewVC logotype

View of /pkg/tests/poMatrix.Rout.save

Parent Directory Parent Directory | Revision Log Revision Log


Revision 473 - (download) (annotate)
Tue Feb 1 14:41:09 2005 UTC (14 years, 6 months ago) by maechler
File size: 3921 byte(s)
more po-tests
R : Copyright 2005, The R Foundation for Statistical Computing
Version 2.1.0 Under development (unstable) (2005-02-01), ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for a HTML browser interface to help.
Type 'q()' to quit R.

> ### Testing positive definite matrices
> 
> library(Matrix)
> 
> stopifnot(c(0,0) == dim(Hilbert(0)))
> 
> h9 <- Hilbert(9)
> str(h9)
Formal class 'poMatrix' [package "Matrix"] with 5 slots
  ..@ uplo         : chr "U"
  ..@ x            : num [1:81] 1.000 0.500 0.333 0.250 0.200 ...
  ..@ Dim          : int [1:2] 9 9
  ..@ rcond        : num(0) 
  ..@ factorization: list()
> all.equal(determinant(h9)$modulus, -96.7369450737858, tol= 1e-15)
[1] TRUE
> stopifnot(0 == length(h9@factorization))# nothing yet
> str(f9 <- as(chol(h9), "trMatrix"))
Formal class 'trMatrix' [package "Matrix"] with 6 slots
  ..@ uplo         : chr "U"
  ..@ diag         : chr "N"
  ..@ x            : num [1:81] 1.000 0.500 0.333 0.250 0.200 ...
  ..@ Dim          : int [1:2] 9 9
  ..@ rcond        : num(0) 
  ..@ factorization: list()
> ## h9 now has factorization
> stopifnot(names(h9@factorization) == "Cholesky")
> rcond(h9)
[1] 9.093796e-13
> rcond(f9)
[1] 9.127195e-07
> str(h9)# has 'rcond' and 'factorization'
Formal class 'poMatrix' [package "Matrix"] with 5 slots
  ..@ uplo         : chr "U"
  ..@ x            : num [1:81] 1.000 0.500 0.333 0.250 0.200 ...
  ..@ Dim          : int [1:2] 9 9
  ..@ rcond        : Named num 9.1e-13
  .. ..- attr(*, "names")= chr "O"
  ..@ factorization:List of 1
  .. ..$ Cholesky:Formal class 'Cholesky' [package "Matrix"] with 6 slots
  .. .. .. ..@ uplo         : chr "U"
  .. .. .. ..@ diag         : chr "N"
  .. .. .. ..@ x            : num [1:81] 1.000 0.500 0.333 0.250 0.200 ...
  .. .. .. ..@ Dim          : int [1:2] 9 9
  .. .. .. ..@ rcond        : num(0) 
  .. .. .. ..@ factorization: list()
> options(digits=4)
> crossprod(f9)# looks the same as
        [,1]   [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]    [,9]
 [1,] 1.0000 0.5000 0.33333 0.25000 0.20000 0.16667 0.14286 0.12500 0.11111
 [2,] 0.5000 0.3333 0.25000 0.20000 0.16667 0.14286 0.12500 0.11111 0.10000
 [3,] 0.3333 0.2500 0.20000 0.16667 0.14286 0.12500 0.11111 0.10000 0.09091
 [4,] 0.2500 0.2000 0.16667 0.14286 0.12500 0.11111 0.10000 0.09091 0.08333
 [5,] 0.2000 0.1667 0.14286 0.12500 0.11111 0.10000 0.09091 0.08333 0.07692
 [6,] 0.1667 0.1429 0.12500 0.11111 0.10000 0.09091 0.08333 0.07692 0.07143
 [7,] 0.1429 0.1250 0.11111 0.10000 0.09091 0.08333 0.07692 0.07143 0.06667
 [8,] 0.1250 0.1111 0.10000 0.09091 0.08333 0.07692 0.07143 0.06667 0.06250
 [9,] 0.1111 0.1000 0.09091 0.08333 0.07692 0.07143 0.06667 0.06250 0.05882
> h9 # but not internally!
        [,1]   [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]    [,9]
 [1,] 1.0000 0.5000 0.33333 0.25000 0.20000 0.16667 0.14286 0.12500 0.11111
 [2,] 0.5000 0.3333 0.25000 0.20000 0.16667 0.14286 0.12500 0.11111 0.10000
 [3,] 0.3333 0.2500 0.20000 0.16667 0.14286 0.12500 0.11111 0.10000 0.09091
 [4,] 0.2500 0.2000 0.16667 0.14286 0.12500 0.11111 0.10000 0.09091 0.08333
 [5,] 0.2000 0.1667 0.14286 0.12500 0.11111 0.10000 0.09091 0.08333 0.07692
 [6,] 0.1667 0.1429 0.12500 0.11111 0.10000 0.09091 0.08333 0.07692 0.07143
 [7,] 0.1429 0.1250 0.11111 0.10000 0.09091 0.08333 0.07692 0.07143 0.06667
 [8,] 0.1250 0.1111 0.10000 0.09091 0.08333 0.07692 0.07143 0.06667 0.06250
 [9,] 0.1111 0.1000 0.09091 0.08333 0.07692 0.07143 0.06667 0.06250 0.05882
> ##   i.e. this is all wrong : all.equal(h9, crossprod(f9))
> 

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge