[R] qr.solve (fwd)

[Torsten Hothorn]

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