Skip to content

lower.tail option in pnorm

6 messages · Ken Knoblauch, Brian Ripley, Ben Bolker +2 more

Original
Ken Knoblauch
# 
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)
[1] -0.6931472 -1.8410216 -0.6931472 -0.1727538 -0.6931472
  [6] -0.1727538 -0.6931472 -1.8410216 -0.6931472 -1.8410216
[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.
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
Brian Ripley
# 
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:

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Ben Bolker
    

    

      
      #
    

  


  

      
Ken Knoblauch <ken.knoblauch <at> inserm.fr> writes:

            

      
      
      
    

        
[snip]

            

      
      
      
    

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

  
  


  


      

    

    

      
      

        


  

        

      Ken Knoblauch
    

    

      
      #
    

  


  

      
Thank you.  That explains it.  I didn't read closely enough.

best,

Ken

Quoting Prof Brian Ripley <ripley at stats.ox.ac.uk>:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Ravi Varadhan
    

    

      
      #
    

  


  

      
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

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      William Dunlap
    

    

      
      #
    

  


  

      
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