Dear All,
I have written the following programs to find a non-singular (10*10)
covariance matrix.
Here is the program:
nitems <- 10
x <- array(rnorm(5*nitems,3,3), c(5,nitems))
sigma <- t(x)%*%x
inverse <- try(solve(sigma), TRUE)
while(inherits(inverse, "try-error"))
{
x <- array(rnorm(5*nitems,3,3), c(5,nitems))
sigma <- t(x)%*%x
inverse <- try(solve(sigma), TRUE)
}
The loop doesn't stop ... This means that no "non-singular"
matrix was found!!!
some thing wrong !!
Thanks a lot for any reply
Bernard
---------------------------------
[[alternative HTML version deleted]]
Don't you have the dimension of x backward? Try:> set.seed(1) > x <- matrix(rnorm(50, 3, 3), 10, 5) > vinv <- solve(crossprod(x)) > vinv[,1] [,2] [,3] [,4] [,5] [1,] 0.019918251 -0.006247646 0.006600209 0.003687249 -0.018670806 [2,] -0.006247646 0.018121025 -0.014815905 -0.005647350 0.003434065 [3,] 0.006600209 -0.014815905 0.023411617 -0.002250342 -0.003258960 [4,] 0.003687249 -0.005647350 -0.002250342 0.025168959 -0.020070844 [5,] -0.018670806 0.003434065 -0.003258960 -0.020070844 0.039593016 If you really have 5 cases and 10 variables, the covariance matrix will have to be singular. Andy> From: Marc Bernard > > Dear All, > > I have written the following programs to find a > non-singular (10*10) covariance matrix. > Here is the program: > > nitems <- 10 > > x <- array(rnorm(5*nitems,3,3), c(5,nitems)) > > sigma <- t(x)%*%x > > inverse <- try(solve(sigma), TRUE) > > > > while(inherits(inverse, "try-error")) > > { > > x <- array(rnorm(5*nitems,3,3), c(5,nitems)) > > sigma <- t(x)%*%x > > inverse <- try(solve(sigma), TRUE) > > } > > > > The loop doesn't stop ... This means that no "non-singular" > matrix was found!!! > > some thing wrong !! > > > > Thanks a lot for any reply > > > > Bernard > > > > > > > > > > --------------------------------- > > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! > http://www.R-project.org/posting-guide.html > >
'x' is a rank 5 matrix from which you want to create a rank 10
crossproduct! Try the following instead:
x <- array(rnorm(100 * nitems, 3, 3), c(100, nitems))
sigma <- crossprod(x)
I hope it helps.
Best,
Dimitris
----
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven
Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/(0)16/336899
Fax: +32/(0)16/337015
Web: http://www.med.kuleuven.be/biostat/
http://www.student.kuleuven.be/~m0390867/dimitris.htm
----- Original Message -----
From: "Marc Bernard" <bernarduse1 at yahoo.fr>
To: <r-help at stat.math.ethz.ch>
Sent: Thursday, October 06, 2005 3:11 PM
Subject: [R] Singular matrix
> Dear All,
>
> I have written the following programs to find a non-singular
> (10*10) covariance matrix.
> Here is the program:
>
> nitems <- 10
>
> x <- array(rnorm(5*nitems,3,3), c(5,nitems))
>
> sigma <- t(x)%*%x
>
> inverse <- try(solve(sigma), TRUE)
>
>
>
> while(inherits(inverse, "try-error"))
>
> {
>
> x <- array(rnorm(5*nitems,3,3), c(5,nitems))
>
> sigma <- t(x)%*%x
>
> inverse <- try(solve(sigma), TRUE)
>
> }
>
>
>
> The loop doesn't stop ... This means that no "non-singular"
matrix
> was found!!!
>
> some thing wrong !!
>
>
>
> Thanks a lot for any reply
>
>
>
> Bernard
>
>
>
>
>
>
>
>
>
> ---------------------------------
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
> http://www.R-project.org/posting-guide.html
>
Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
>>>>> "Marc" == Marc Bernard <bernarduse1 at yahoo.fr> >>>>> on Thu, 6 Oct 2005 15:11:09 +0200 (CEST) writes:........... Marc> Here is the program: Marc> nitems <- 10 Marc> x <- array(rnorm(5*nitems,3,3), c(5,nitems)) Marc> sigma <- t(x)%*%x Marc> inverse <- try(solve(sigma), TRUE) ............. Just a side remark on your code above: You should learn about and use crossprod(x) rather than t(x) %*% x because of a 1) more efficient implementation 2) more accurate implementation The exact details depends on the (accelerated/optimized or not) version BLAS/Lapack your version of R is using and also on the kind of matrices. Further note, that for crossprod(t(X)) == X %*% t(X) LAPACK also provides a direc version to which the 'Matrix' package interfaces via function tcrossprod() which in particular also works *fast* for some of the sparse matrix classes. Martin Maechler, ETH Zurich