[R] qr.solve (fwd)

[Thomas Lumley]

On Tue, 14 Mar 2000, Torsten Hothorn wrote:

> 
> Two friend reported me a problem, which I can't solve:
> 
> (I run R-1.0.0, Debian Linux) 
> 
> They hava a function "corr.matrix" (see end of mail), and when they
> create a 173x173 matrix with this function
> 
> V <- corr.matrix(0.3, n=173)
> V1 <- qr.solve(V)
> 
> reports:
> 
> Error in qr(a, tol = tol) : NA/NaN/Inf in foreign function call (arg 1)
> 
> For n < 173, qr.solve returns the correct result. 


Well, that's because there _are_ NaNs in V (two of them). There are also
some incorrect zeros due to floating point overflow.

In computing V you calculate gamma(172), which is a really really big
number (about 10^309) and it overflows to Inf. As a result, V[1,172]=0 and
V[1,173]=NaN.  

Incidentally, for n=172 you still don't get the right answer.  The
overflow occurs only in the denominator and the [1,172] element is 0
rather than about 0.05.

Use lgamma instead of gamma and exponentiate the result. 

	-thomas




> 
> Torsten
> 
> ________________________________________________________________________
> 
> corr.matrix <- function(d, n=100)
> {
>   mat <- NULL
>   for(i in 1:n)
>   {
>   mat <- c(mat,((1-i):(n-i)))
>   }
>   mat <- abs(mat)
>   mat <- corr(mat,d)
>   mat <- matrix(mat, ncol=n)
>   mat
> }
> 
> corr <- function (x,d)
> {
>   rho <- gamma(x+d)*gamma(1-d)
>   rho <- rho/(gamma(x-d+1)*gamma(d))
>   rho    
> }
> 
> tr <- function(M)
> {
>   sum(diag(M))
> }
> 
> eff <- function(d, n=100, tol=1e-7)
> {
>   x <- 1:n
>   X <- cbind(1,x)
>   V <- corr.matrix(d, n=n)
>   V1 <- qr.solve(V, tol=tol)
>   GLS <- (t(X)%*%V1%*%X)
>   GLS <- solve(GLS) %*%t(X)%*%X
>   OLS <- solve(t(X)%*%X)%*%t(X)%*%V%*%X
>   eff <- tr(GLS)/tr(OLS)
>   eff
> }
> 
> eff.all <- function(n=100, int=0.01, tol=1e-7)
> {
>   d <- -0.49
>   x <- NULL
>   y <- NULL
>   while(d<0.5)
>   {
>     x <- c(x,d)
>     y <- c(y,eff(d,n=n))
>     d <- d + int
>   }   
>   plot(x,y, type="l", xlab="d", ylab="eff (T, d)", main="Relative efficiency of OLS")
>   d <- x
>   effizienz <- y
>   result <- list(d,effizienz)
>   names(result) <- c("d", "effizienz")
>   return(result)
> }
> 

Thomas Lumley
Assistant Professor, Biostatistics
University of Washington, Seattle

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._