Hi, I would have thought that these two constructions would produce the same result but they do not. Resp <- rbinom(10, 1, 0.5) Stim <- rep(0:1, 5) mm <- model.matrix(~ Stim) Xb <- mm %*% c(0, 1) ifelse(Resp, log(pnorm(Xb)), log(1 - pnorm(Xb))) pnorm(as.vector(Xb), lower.tail = Resp, log.p = TRUE)> ifelse(Resp, log(pnorm(Xb)), log(1 - pnorm(Xb)))[1] -0.6931472 -1.8410216 -0.6931472 -0.1727538 -0.6931472 [6] -0.1727538 -0.6931472 -1.8410216 -0.6931472 -1.8410216> pnorm(as.vector(Xb), lower.tail = Resp, log.p = TRUE)[1] -0.6931472 -1.8410216 -0.6931472 -1.8410216 -0.6931472 [6] -1.8410216 -0.6931472 -1.8410216 -0.6931472 -1.8410216 If I have missed something obvious, I would be grateful to have it pointed out.> sessionInfo()R version 2.10.1 beta (2009-12-04 r50668) i386-apple-darwin9.8.0 locale: [1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] stats graphics grDevices utils datasets methods [7] base loaded via a namespace (and not attached): [1] tools_2.10.1 Thanks for any enlightenment. best, Ken -- Ken Knoblauch Inserm U846 Stem-cell and Brain Research Institute Department of Integrative Neurosciences 18 avenue du Doyen L?pine 69500 Bron France tel: +33 (0)4 72 91 34 77 fax: +33 (0)4 72 91 34 61 portable: +33 (0)6 84 10 64 10 http://www.sbri.fr/members/kenneth-knoblauch.html
>From the help page:pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE) x,q: vector of quantiles. lower.tail: logical; if TRUE (default), probabilities are P[X <= x], otherwise, P[X > x]. Note that lower.tail is not said to be a vector, and the first value is taken (what it is is random in your example). On Tue, 8 Dec 2009, Ken Knoblauch wrote:> Hi, > > I would have thought that these two constructions would > produce the same result but they do not. > > Resp <- rbinom(10, 1, 0.5) > Stim <- rep(0:1, 5) > mm <- model.matrix(~ Stim) > Xb <- mm %*% c(0, 1) > ifelse(Resp, log(pnorm(Xb)), log(1 - pnorm(Xb))) > pnorm(as.vector(Xb), lower.tail = Resp, log.p = TRUE) > >> ifelse(Resp, log(pnorm(Xb)), log(1 - pnorm(Xb))) > [1] -0.6931472 -1.8410216 -0.6931472 -0.1727538 -0.6931472 > [6] -0.1727538 -0.6931472 -1.8410216 -0.6931472 -1.8410216 >> pnorm(as.vector(Xb), lower.tail = Resp, log.p = TRUE) > [1] -0.6931472 -1.8410216 -0.6931472 -1.8410216 -0.6931472 > [6] -1.8410216 -0.6931472 -1.8410216 -0.6931472 -1.8410216 > > If I have missed something obvious, I would be grateful > to have it pointed out. > >> sessionInfo() > R version 2.10.1 beta (2009-12-04 r50668) > i386-apple-darwin9.8.0 > > locale: > [1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8 > > attached base packages: > [1] stats graphics grDevices utils datasets methods > [7] base > > loaded via a namespace (and not attached): > [1] tools_2.10.1 > > Thanks for any enlightenment. > > best, > > Ken > > -- > Ken Knoblauch > Inserm U846 > Stem-cell and Brain Research Institute > Department of Integrative Neurosciences > 18 avenue du Doyen L?pine > 69500 Bron > France > tel: +33 (0)4 72 91 34 77 > fax: +33 (0)4 72 91 34 61 > portable: +33 (0)6 84 10 64 10 > http://www.sbri.fr/members/kenneth-knoblauch.html > > ______________________________________________ > R-devel at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel-- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595
Ken Knoblauch <ken.knoblauch <at> inserm.fr> writes:> > Hi, > > I would have thought that these two constructions would > produce the same result but they do not. > > Resp <- rbinom(10, 1, 0.5) > Stim <- rep(0:1, 5) > mm <- model.matrix(~ Stim) > Xb <- mm %*% c(0, 1) > ifelse(Resp, log(pnorm(Xb)), log(1 - pnorm(Xb))) > pnorm(as.vector(Xb), lower.tail = Resp, log.p = TRUE) >[snip]> If I have missed something obvious, I would be grateful > to have it pointed out. >lower.tail is not vectorized. All elements but the first are ignored. This seems fairly obvious to me from reading ?pnorm (e.g. "mean" is described a "vector of means", "sd" is described as "vector of standard deviations", but lower.tail is described as "logical", but 'obvious' is certainly in the eye of the beholder. You can do what you want with mapply(pnorm,q=c(-1,1),lower.tail=c(0,1)) .
They will not be the same. The problem is that the `lower.tail' argument is not vectorized. Therefore, it is always set equal to the first element of `Resp', which in your example is FALSE. If you want to obtain same results, this will do the trick: ans1 <- ifelse(Resp, log(pnorm(Xb)), log(1 - pnorm(Xb))) ans2 <- apply(cbind(Xb, Resp), 1, function(x) pnorm(x[1], lower.tail=x[2], log.p=TRUE)) all.equal(ans1, as.numeric(ans2)) Ravi. ____________________________________________________________________ Ravi Varadhan, Ph.D. Assistant Professor, Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins University Ph. (410) 502-2619 email: rvaradhan at jhmi.edu ----- Original Message ----- From: Ken Knoblauch <ken.knoblauch at inserm.fr> Date: Tuesday, December 8, 2009 10:37 am Subject: [Rd] lower.tail option in pnorm To: r-devel at r-project.org> Hi, > > I would have thought that these two constructions would > produce the same result but they do not. > > Resp <- rbinom(10, 1, 0.5) > Stim <- rep(0:1, 5) > mm <- model.matrix(~ Stim) > Xb <- mm %*% c(0, 1) > ifelse(Resp, log(pnorm(Xb)), log(1 - pnorm(Xb))) > pnorm(as.vector(Xb), lower.tail = Resp, log.p = TRUE) > > > ifelse(Resp, log(pnorm(Xb)), log(1 - pnorm(Xb))) > [1] -0.6931472 -1.8410216 -0.6931472 -0.1727538 -0.6931472 > [6] -0.1727538 -0.6931472 -1.8410216 -0.6931472 -1.8410216 > > pnorm(as.vector(Xb), lower.tail = Resp, log.p = TRUE) > [1] -0.6931472 -1.8410216 -0.6931472 -1.8410216 -0.6931472 > [6] -1.8410216 -0.6931472 -1.8410216 -0.6931472 -1.8410216 > > If I have missed something obvious, I would be grateful > to have it pointed out. > > > sessionInfo() > R version 2.10.1 beta (2009-12-04 r50668) > i386-apple-darwin9.8.0 > > locale: > [1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8 > > attached base packages: > [1] stats graphics grDevices utils datasets methods > [7] base > > loaded via a namespace (and not attached): > [1] tools_2.10.1 > > Thanks for any enlightenment. > > best, > > Ken > > -- > Ken Knoblauch > Inserm U846 > Stem-cell and Brain Research Institute > Department of Integrative Neurosciences > 18 avenue du Doyen L?pine > 69500 Bron > France > tel: +33 (0)4 72 91 34 77 > fax: +33 (0)4 72 91 34 61 > portable: +33 (0)6 84 10 64 10 > > > ______________________________________________ > R-devel at r-project.org mailing list >
A more direct way to reproduce this is > pbinom(q=0:4, size=4, prob=.25, lower.tail=FALSE) [1] 0.68359375 0.26171875 0.05078125 0.00390625 0.00000000 > pbinom(q=0:4, size=4, prob=.25, lower.tail=c(FALSE,TRUE,TRUE,TRUE)) [1] 0.68359375 0.26171875 0.05078125 0.00390625 0.00000000 > pbinom(q=0:4, size=4, prob=.25, lower.tail=c(TRUE,TRUE,TRUE,TRUE)) [1] 0.31640625 0.73828125 0.94921875 0.99609375 1.00000000 which shows that the lower.tail argument is not not vectorized. While this fact may be documented, it would be nice if functions that expected a scalar argument would complain if they got a nonscalar argument. Some do complain, as in > seq(1, 1:4) [1] 1 Warning message: In from:to : numerical expression has 4 elements: only the first used but lots silently take the first element of a vector when its length is more than 1, as above, and use the default value if the length is 0, as in > pbinom(q=0:4, size=4, prob=.25, lower.tail=logical()) [1] 0.31640625 0.73828125 0.94921875 0.99609375 1.00000000 Should the C functions Rf_asReal, Rf_asInteger, etc., warn if the argument has length more than 1? If so, they probably need another argument, a string, to put into the warning message, so the user knows which argument is causing the complaint. Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com> -----Original Message----- > From: r-devel-bounces at r-project.org > [mailto:r-devel-bounces at r-project.org] On Behalf Of Ken Knoblauch > Sent: Tuesday, December 08, 2009 7:36 AM > To: r-devel at r-project.org > Subject: [Rd] lower.tail option in pnorm > > Hi, > > I would have thought that these two constructions would > produce the same result but they do not. > > Resp <- rbinom(10, 1, 0.5) > Stim <- rep(0:1, 5) > mm <- model.matrix(~ Stim) > Xb <- mm %*% c(0, 1) > ifelse(Resp, log(pnorm(Xb)), log(1 - pnorm(Xb))) > pnorm(as.vector(Xb), lower.tail = Resp, log.p = TRUE) > > > ifelse(Resp, log(pnorm(Xb)), log(1 - pnorm(Xb))) > [1] -0.6931472 -1.8410216 -0.6931472 -0.1727538 -0.6931472 > [6] -0.1727538 -0.6931472 -1.8410216 -0.6931472 -1.8410216 > > pnorm(as.vector(Xb), lower.tail = Resp, log.p = TRUE) > [1] -0.6931472 -1.8410216 -0.6931472 -1.8410216 -0.6931472 > [6] -1.8410216 -0.6931472 -1.8410216 -0.6931472 -1.8410216 > > If I have missed something obvious, I would be grateful > to have it pointed out. > > > sessionInfo() > R version 2.10.1 beta (2009-12-04 r50668) > i386-apple-darwin9.8.0 > > locale: > [1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8 > > attached base packages: > [1] stats graphics grDevices utils datasets methods > [7] base > > loaded via a namespace (and not attached): > [1] tools_2.10.1 > > Thanks for any enlightenment. > > best, > > Ken > > -- > Ken Knoblauch > Inserm U846 > Stem-cell and Brain Research Institute > Department of Integrative Neurosciences > 18 avenue du Doyen L?pine > 69500 Bron > France > tel: +33 (0)4 72 91 34 77 > fax: +33 (0)4 72 91 34 61 > portable: +33 (0)6 84 10 64 10 > http://www.sbri.fr/members/kenneth-knoblauch.html > > ______________________________________________ > R-devel at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel >