thr3ads.net - R devel - [Rd] lower.tail option in pnorm [Dec 2009]

If this information is useful, please help other people find it:
Share via:

Ken Knoblauch

2009-Dec-08 15:36 UTC

[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

Prof Brian Ripley

2009-Dec-08 16:20 UTC

head link

[Rd] lower.tail option in pnorm

>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

Ben Bolker

2009-Dec-08 16:23 UTC

head link

[Rd] lower.tail option in pnorm

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))
.

Ravi Varadhan

2009-Dec-08 16:32 UTC

head link

[Rd] lower.tail option in pnorm

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
>

William Dunlap

2009-Dec-08 16:46 UTC

head link

[Rd] lower.tail option in pnorm

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
>

Reasonably Related Threads

Search for more possibly parallel threads

R devel - Dec 2009 - lower.tail option in pnorm

[Rd] lower.tail option in pnorm

[Rd] lower.tail option in pnorm

[Rd] lower.tail option in pnorm

[Rd] lower.tail option in pnorm

[Rd] lower.tail option in pnorm

Reasonably Related Threads