[R] quantreg::rq.fit.hogg crashing at random

[Madison Lincoln]

Dear all,

I am using the "rq.fit.hogg" function from the "quantreg" package. I have
two problems with it.

First (less importantly), it gives an error at its default values with
error message "Error in if (n2 != length(r)) stop("R and r of incompatible
dimension") : argument is of length zero". I solve this by commenting four
lines in the code. I.e. I define a new function "rq.fit.hogg2" that is the
same as "rq.fit.hogg" but with four lines commented. You will see this in
my code at the end of this post. I understand why this solution works, so
it is not really a problem right now.

Second (importantly), "rq.fit.hogg" crashes frequently. The message I get
is "R for Windows GUI front-end has stopped working. A problem cause the
program to stop working correctly. Windows will close the program and
notify you if a solution is available."

I don't know how to provide a reproducible example of the crash because the
function seems to crash at random. That is, I generate some data (with a
fixed seed so that I can replicate it exactly), input into the function,
and sometimes it crashes but sometimes not. If I create a loop over
different seeds and run the function once in each iteration, sometimes it
crashes on the first iteration while sometimes it goes fine for some 20
iterations and crashes only then.

I am including the code with the loop that should eventually produce a
crash. If it does not, try running it again; that worked every time for me.
I am also including session info and locale info.

Please help me. Thank you!

Kind regards,


*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

Here are my session details etc.:

> Sys.getlocale()
[1] "LC_COLLATE=English_United States.1252;LC_CTYPE=English_United
States.1252;LC_MONETARY=English_United
States.1252;LC_NUMERIC=C;LC_TIME=English_United
States.1252"

> sessionInfo()
R version 3.4.0 (2017-04-21)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 14393)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252
[2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] quantreg_5.33 SparseM_1.77  MASS_7.3-47

loaded via a namespace (and not attached):
[1] compiler_3.4.0     Matrix_1.2-10      MatrixModels_0.4-1 grid_3.4.0

[5] lattice_0.20-35

*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

Here is the code that crashes:

#-----Here I redefine the function as described in my first paragraph

rq.fit.hogg2 = function (x, y, taus = c(0.1, 0.3, 0.5), weights = c(0.7,
0.2, 0.1), R = NULL, r = NULL, beta = 0.99995, eps = 1e-06)
{
    n <- length(y)
    n2 <- nrow(R)
    m <- length(taus)
    p <- ncol(x) + m
    if (n != nrow(x))
        stop("x and y don't match n")
    if (m != length(weights))
        stop("taus and weights differ in length")
    if (any(taus < eps) || any(taus > 1 - eps))
        stop("taus outside (0,1)")
    W <- diag(weights)
    if (m == 1)
        W <- weights
    x <- as.matrix(x)
    X <- cbind(kronecker(W, rep(1, n)), kronecker(weights, x))
    y <- kronecker(weights, y)
    rhs <- c(weights * (1 - taus) * n, sum(weights * (1 - taus)) *
        apply(x, 2, sum))
#    if (n2 != length(r))
#        stop("R and r of incompatible dimension")
#    if (ncol(R) != p)
#        stop("R and X of incompatible dimension")
    d <- rep(1, m * n)
    u <- rep(1, m * n)
    if (length(r)) {
        wn1 <- rep(0, 10 * m * n)
        wn1[1:(m * n)] <- 0.5
        wn2 <- rep(0, 6 * n2)
        wn2[1:n2] <- 1
        z <- .Fortran("rqfnc", as.integer(m * n), as.integer(n2),
            as.integer(p), a1 = as.double(t(as.matrix(X))), c1 =
as.double(-y),
            a2 = as.double(t(as.matrix(R))), c2 = as.double(-r),
            rhs = as.double(rhs), d1 = double(m * n), d2 = double(n2),
            as.double(u), beta = as.double(beta), eps = as.double(eps),
            wn1 = as.double(wn1), wn2 = as.double(wn2), wp = double((p +
                3) * p), it.count = integer(3), info = integer(1))
    }
    else {
        wn <- rep(0, 10 * m * n)
        wn[1:(m * n)] <- 0.5
        z <- .Fortran("rqfnb", as.integer(m * n), as.integer(p),
            a = as.double(t(as.matrix(X))), c = as.double(-y),
            rhs = as.double(rhs), d = as.double(d), as.double(u),
            beta = as.double(beta), eps = as.double(eps), wn =
as.double(wn),
            wp = double((p + 3) * p), it.count = integer(2),
            info = integer(1))
    }
    if (z$info != 0)
        warning(paste("Info = ", z$info, "in stepy: singular design:
iterations ",
            z$it.count[1]))
    coefficients <- -z$wp[1:p]
    if (any(is.na(coefficients)))
        stop("NA coefs:  infeasible problem?")
    list(coefficients = coefficients, nit = z$it.count, flag = z$info)
}

#----- Here is the main program that crashes

library(MASS)
library(quantreg)

M=1e3; n=1e3; p=15; type=8; K=10 # no. of replications, sample size,
dimension of beta, type of quantile estimator, number of quantiles
distributions=c("norm","t1","t2","t3","t5","logistic","
exponential","Weibull")
method.wqr="fn" # interior point method; if "fn" then roughly matches with
method.cqr="ip"

# Create the covariance matrix of X
Sigma=matrix(NA,p,p); for(i in 1:p) for(j in 1:p) Sigma[i,j]=0.5^(abs(i-j))
# Generate X (common across all simulations)
set.seed(0); X=mvrnorm(n=n,mu=rep(0,p),Sigma=Sigma)

Binvlist=list()
for(k in 1:K){
 tau=cumsum(rep(1/(k+1),k))
 Ai=matrix(rep(tau,k),nrow=k,ncol=k,byrow=TRUE)
 Aj=matrix(rep(tau,k),nrow=k,ncol=k,byrow=FALSE)
 Amin=pmin(Ai,Aj) # Amin=Ai; Amin[Ai>Aj]=Aj[Ai>Aj]
 Ax=tau %*% t(tau)
 B=Amin-Ax
 Binvlist[[k]]=solve(B)
}

for(m in 1:M){

mse_wqr_list=mse_cqr_list=list()
j=0; for(distr in distributions){

 j=j+1
 #print(distributions[col])
 mse_wqr_list[[j]]=mse_cqr_list[[j]]=matrix(NA,M,K)
 set.seed(m); beta=rnorm(p)
 set.seed(m+1000000000)
 if(distr=="norm"       ) eps=rnorm(n)
 if(substr(x=distr,start=1,stop=1)=="t"){
df=as.numeric(substr(x=distr,start=2,stop=nchar(distr)));
eps=rt(n,df=df) }
 if(distr=="logistic"   ) eps=rlogis(n)
 if(distr=="exponential") eps=rexp(n)
 if(distr=="Weibull"    ) eps=rweibull(n,shape=1.5)
 y=X%*%beta+eps
 model=lm(y~X); res=model$res; d=density(res)

 mse_wqr=mse_cqr=rep(NA,K)
 for(k in 1:K){

  #----- Find estimated optimal weights
  tau=cumsum(rep(1/(k+1),k))
  Binv=Binvlist[[k]]
   q = quantile(res,probs=tau,type=type) # Hyndman recommends `type=8` for
the `quantile` function
   tmp=as.numeric(c(d$x,q)); ranks=rank(tmp);
below=tail(ranks,length(q))-c(1:length(q));
above=below+1
   v = d$y[below] * (q-d$x[below])/(d$x[above]-d$x[below]) + d$y[above] *
(d$x[above]-q)/(d$x[above]-d$x[below])
   if(k==1) V = matrix(v,1,1) else V = diag(v)
  w_CQR = Binv %*% v; w_CQR = w_CQR / sum(w_CQR)

  #----- Use the estimated optimal weights to actually estimate beta (with
EQR and CQR) and evaluate how well it does
  print(paste(m, distr, k)); readline()
  ################ THE NEXT LINE CRASHES AT RANDOM:
  fit_cqr=try( rq.fit.hogg2(x=cbind(X),y=y,taus=tau,weights=c(w_CQR)) )

 } # for(k in 1:K)

} # for(distr in distributions)

} # for(m in 1:M)

	[[alternative HTML version deleted]]