# [R] qr.solve (fwd)

Torsten Hothorn hothorn at ci.tuwien.ac.at
Tue Mar 14 17:22:19 CET 2000

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