[R-SIG-Mac] dblepr for 64B R

[Roger Koenker]

Hi David,

Thanks, but this picture is exactly what I was getting in 64B.....   
and it is entirely bogus.
Fortunately, this morning I discovered that ddot in fortran needs to  
be declared double
in 64B, but apparently doesn't in 32B... and this resolves the  
problem.  I've submitted
an updated version of quantreg to CRAN that fixes the problem.   The  
dblepr was just
an intermediate debugging issue, and was entirely my stupidity.  It  
works exactly
as it should once variables are declared double.  Why ddot behaves  
differently in 64B
vs 32B is still mysterious to me, but so are many aspects of ordinary  
life, so I have
resolved not to question further.


Roger Koenker
rkoenker at illinois.edu




On Jun 17, 2010, at 8:00 PM, David Winsemius wrote:

>
> On Jun 16, 2010, at 5:24 PM, Roger Koenker wrote:
>
>> I'm trying to understand why 64B R produces an insane answer/plot
>> to the first
>>
>> 	example(crq)
>>
>> in my quantreg package, whereas 32B R produces a sane answer
>> with the same data.  This involves some fortran and in my usual
>> primitive fashion I've been trying to debug with fortran printing
>> using the tried and true call dblepr.  However, in my test function
>> call dblepr prints  values like
>> z
>> [1] 1.644287e-313
>>
>> even though the function returns the sane, correct values.  Is this
>> a known issue, and if so is there an alternative strategy for
>> printing from fortran in 64b R?
>>
>> TIA for any enlightenment.
>> Roger
>>
>> PS.  In the hope that someone will tell me that all would be well
>> if I only followed directions and upgraded R,  I'll sheepishly admit:
>>
>>> sessionInfo()
>> R version 2.11.0 Under development (unstable) (2010-02-18 r51149)
>> x86_64-apple-darwin9.8.0
>
> With a freshly installed 2.11.1 just a couple of days ago, and I  
> only run the 64-bit version, and I do not see any such result. There  
> is a line on the first pot that "shoots up". I attached  
> it.<First.quantreg.example.plot.pdf>
> Here is the console output I get:
>
> > example(crq)
>
> crq> # An artificial Powell example
> crq> set.seed(2345)
>
> crq> x <- sqrt(rnorm(100)^2)
>
> crq> y <-  -0.5 + x +(.25 + .25*x)*rnorm(100)
>
> crq> plot(x,y, type="n")
> Hit <Return> to see next plot:
>
> crq> s <- (y > 0)
>
> crq> points(x[s],y[s],cex=.9,pch=16)
>
> crq> points(x[!s],y[!s],cex=.9,pch=1)
>
> crq> yLatent <- y
>
> crq> y <- pmax(0,y)
>
> crq> yc <- rep(0,100)
>
> crq> for(tau in (1:4)/5){
> crq+         f <- crq(Curv(y,yc) ~ x, tau = tau, method = "Pow")
> crq+         xs <- sort(x)
> crq+         lines(xs,pmax(0,cbind(1,xs)%*%f$coef),col="red")
> crq+         abline(rq(y ~ x, tau = tau), col="blue")
> crq+         abline(rq(yLatent ~ x, tau = tau), col="green")
> crq+         }
> Loading required package: survival
> Loading required package: splines
>
> Attaching package: 'survival'
>
> The following object(s) are masked from 'package:quantreg':
>
>    untangle.specials
>
>
> crq> legend(.15,2.5,c("Naive QR","Censored QR","Omniscient QR"),
> crq+         lty=rep(1,3),col=c("blue","red","green"))
>
> crq> data(uis)
>
> crq> #estimate the Peng and Huang model using log(TIME) AFT  
> specification
> crq> fit <- crq(Surv(log(TIME), CENSOR) ~  ND1 + ND2 + IV3 +
> crq+                TREAT + FRAC + RACE + AGE * SITE, method =  
> "Por", data = uis)
>
> crq> Sfit <- summary(fit,1:19/20)
> bootstrap roughly  10  percent complete
> bootstrap roughly  20  percent complete
> bootstrap roughly  30  percent complete
> bootstrap roughly  40  percent complete
> bootstrap roughly  50  percent complete
> bootstrap roughly  60  percent complete
> bootstrap roughly  70  percent complete
> bootstrap roughly  80  percent complete
> bootstrap roughly  90  percent complete
> bootstrap roughly  100  percent complete
>
> crq> PHit <- coxph(Surv(TIME, CENSOR) ~  ND1 + ND2 + IV3 +
> crq+                TREAT + FRAC + RACE + AGE * SITE, data = uis)
>
> crq> plot(Sfit, CoxPHit = PHit)
> Hit <Return> to see next plot:
>
> crq> formula <-  ~  ND1 + ND2 + IV3 + TREAT + FRAC + RACE + AGE *  
> SITE -1
>
> crq> X <- data.frame(model.matrix(formula,data=uis))
>
> crq> newd <- as.list(apply(X,2,median))
>
> crq> pred <- predict(fit, newdata=newd, type = "stepfun")
>
> crq> plot(pred,do.points=FALSE,xlab = expression(tau), ylab =  
> expression(Q(tau)),main= "Quantiles at Median Covariate Values")
> Hit <Return> to see next plot:
>
> crq> plot(rearrange(pred),add=TRUE,do.points=FALSE,col.vert ="red",  
> col.hor="red")
>
> crq> legend(.15,7,c("Raw","Rearranged"),lty =  
> rep(1,2),col=c("black","red"))
> Warning messages:
> 1: In crq.fit.pow(X, y, cen, tau = taus, weights, left = left, ...) :
>  Solution may be nonunique
> 2: In crq.fit.pow(X, y, cen, tau = taus, weights, left = left, ...) :
>  Solution may be nonunique
> 3: In crq.fit.pow(X, y, cen, tau = taus, weights, left = left, ...) :
>  Solution may be nonunique
> 4: In crq.fit.pow(X, y, cen, tau = taus, weights, left = left, ...) :
>  Solution may be nonunique
> >
>
> None of the plots have weird limits. If you tell me what objects to  
> inspect I would be happy to deliver more details. I do not see that  
> dblepr() is in the NAMESPACE.
>
> -- 
> David.
>
>
>>
>>
>> PPS.  The fact that I can run 32B R on the same machine by simply
>> typing R --arch i386  is one of the great advances of modern science
>> in my opinion, or at least I thought so, until I discovered this  
>> discrepancy
>> in what was being computed by the two versions.
>>
>>
>> url:    www.econ.uiuc.edu/~roger            Roger Koenker
>> email    rkoenker at uiuc.edu            Department of Economics
>> vox:     217-333-4558                University of Illinois
>> fax:       217-244-6678                Urbana, IL 61801
>>
>> _______________________________________________
>> R-SIG-Mac mailing list
>> R-SIG-Mac at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-mac
>
> David Winsemius, MD
> West Hartford, CT
>