Dear all,
many thanks to Jon & Ravi for their help on this, and apologies if
r-help would have been a more appropriate forum.
On Tue, 2012-12-11 at 15:43 +0000, Jon Clayden wrote:
> Strategy 1: Some code like this:
> if (det(X) < epsilon) {
> warning("Near singular matrix")
> return(NULL)
> }
> return(solve(X))
>
> This solution is probably the easiest one to take, but to match
> solve.default, the test should be
>
> if (rcond(X) < .Machine$double.eps)
>
> Catching that case should avoid the error. I hope this helps.
Yes, that works well.
On Tue, 2012-12-11 at 15:46 +0000, Ravi Varadhan wrote:
> In any case, you might want to try MASS::ginv instead of solve(), if
> you expect ill-conditioning. Here is one possible solution:
>
> f <- function(X) {
> invX <- tryCatch(ginv(X, tol=.Machine$double.eps),
> error=function(e) {
> warning(e)
> error.flag <- TRUE}) # you should avoid the global
> assignment `<<-'
> if (error.flag) return(NULL)
> return(invX)
> }
I'll try this if any problems do emerge with the above solution. One
small point: I used <<- because there seemed to be no other way of
returning from the function f if an error had been thrown by ginv() (or
solve()):
> flag <- FALSE
> tryCatch(stop(), error=function(e) {flag <- TRUE})
> flag
[1] FALSE
> tryCatch(stop(), error=function(e) {flag <<- TRUE})
> flag
[1] TRUE
I feel sure that there must be a more elegant way of doing this!
David.
--
David C Sterratt, Research Fellow http://homepages.inf.ed.ac.uk/sterratt
Institute for Adaptive and Neural Computation tel: +44 131 651 1739
School of Informatics, University of Edinburgh fax: +44 131 650 6899
Informatics Forum, 10 Crichton Street, Edinburgh EH8 9AB, Scotland
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
NEW BOOK: Principles of Computational Modelling in Neuroscience
Sterratt, Graham, Gillies & Willshaw (CUP, 2011).
http://www.compneuroprinciples.org
The University of Edinburgh is a charitable body, registered in
Scotland, with registration number SC005336.
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
