On 11-02-2013, at 15:00, "Camarda, Carlo Giovanni" <cama...@demogr.mpg.de>
wrote:
> Dear Uwe,
>
> thanks for the response.
> I knew my matrix was almost singular, but is there any way to implement the
> algorithm is solve(base) with a spam-matrix?
>
> As workaround, I might add a ridge penalty:
>
> ApP <- A+10^-1*diag.spam(nrow(A))
> solve(ApP)
>
> but this would completely jeopardize other parts of my model.
I've inspected your matrix with package Matrix.
It doesn't seem catastrophically ill-conditioned.
It's not positive definite.
It seems that spam is issuing a misleading error message.
Try this
ent <- c(2312.12324929972,-2000,1000,-2000,1000,-2000,
5031.91011235955,-2000,-2000,1000,-1,1000,-2000,
2049.8595505618,-2000,1000,-1,-2000,5036.89635854342,
-2000,1000,-2000,1,-2000,-2000,8119.66292134831,-2000,
-2000,-2000,1000,-2000,5058.83426966292,-2000,-1,1000,
-2000,2051.85434173669,-2000,1000,1,1000,-2000,-2000,
5043.87640449438,-2000,1,1000,-2000,1000,-2000,
2110.68820224719,-1,1,0,-1,1,-1,1)
rr <- as.integer(c(1,6,12,18,24,29,35,41,47,52,55,57,59))
cc <- as.integer(c(1,2,3,4,7,1,2,3,5,8,10,1,2,3,6,9,11,1,
4,5,6,7,10,2,4,5,6,8,3,4,5,6,9,12,1,4,
7,8,9,11,2,5,7,8,9,12,3,6,7,8,9,2,4,10,3,7,6,8))
di <- as.integer(c(12,12))
library(Matrix)
A <- sparseMatrix(x=ent,p=rr-1,j=cc)
A
condest(A)
det(A)
solve(A)
kappa(A)
# not pd
# chol(A)
isSymmetric(A)
# chol(A) # ==> not positive definite
chol(A+diag(1e-1,nrow(A))) # this works but isn't what you want
A.m <- as.matrix(A)
chol(A.m,pivot=TRUE) # without pivot: not positive definite
# has rank 9 and warning message
B <- as.matrix(A)
condest(B)
rcond(B)
det(A)
solve(B)
all.equal(as.matrix(A %*% solve(B)), diag(12), check.attributes=FALSE)
all.equal(as.matrix(B %*% solve(A)), diag(12), check.attributes=FALSE)
Berend
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
