# Some functions to debug # John Fox 11 February 05 lregIWLS <- function(X, y, n=rep(1,length(y)), maxIter=10, tol=1E-6){ # bugged! # X is the model matrix # y is the response vector of observed proportion # n is the vector of binomial counts # maxIter is the maximum number of iterations # tol is a convergence criterion X <- cbind(1, X) # add constant b <- bLast <- rep(0, ncol(X)) # initialize it <- 1 # iteration index while (it <= maxIter){ if (max(abs(b - bLast)/(abs(bLast) + 0.01*tol)) < tol) break eta <- X %*% b mu <- 1/(1 + exp(-eta)) nu <- as.vector(mu*(1 - mu)) w <- n*nu z <- eta + (y - mu)/nu b <- lsfit(X, z, w, intercept=FALSE)\$coef bLast <- b it <- it + 1 # increment index } if (it > maxIter) warning('maximum iterations exceeded') Vb <- solve(t(X) %*% diag(w) %*% X) list(coefficients=b, var=Vb, iterations=it) } rowEchelonForm <- function(A){ # bugged! n <- nrow(A) m <- ncol(A) for (i in 1:min(c(m, n))){ currentColumn <- A[,i] currentColumn[1:n < i] <- 0 which <- which.max(abs(currentColumn)) # find maximum pivot in current # column at or below current row pivot <- A[which, i] if (abs(pivot) == 0) next # check for 0 pivot if (which > i) A[c(i, which),] <- A[c(which, i),] # exchange rows A[i,] <- A[i,]/pivot # pivot row <- A[i,] A <- A - outer(A[,i], row) # sweep A[i,] <- row # restore current row } for (i in 1:n) if (max(abs(A[i,1:m])) == 0) A[c(i,n),] <- A[c(n,i),] # 0 rows to bottom A } runningMedian <- function(x, length=3){ # bugged! # x: a numeric vector # length: the number of values for each running median, defaults to 3 n <- length(x) X <- matrix(x, n, length) for (i in 1:length) X[1:(n - i + 1), i] <- x[-(1:(i - 1))] apply(X, 1, median)[1:(n - length + 1)] }