R help - Mar 2000 - qr.solve (fwd)

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. 

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)
}
-------------- next part --------------
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)
}

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