qr.solve (fwd)
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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._