s-news
[Top] [All Lists]

"Bug in qwilcox()?" revisited

To: s-news-list <s-news@lists.biostat.wustl.edu>
Subject: "Bug in qwilcox()?" revisited
From: Gerrit Eichner <Gerrit.Eichner@math.uni-giessen.de>
Date: Tue, 24 Jun 2003 13:57:52 +0200 (CET)
Dear list,

in the meantime I've found out that the strange behaviour of qwilcox() for
m = n = 36 has its origin in the third line of the following segment of
code as it is used in qwilcox():

 x <- ((m * (m + 1))/2):((m * (m + 2 * n + 1))/2)
 cumprob <- cumsum(dwilcox(x, m, n))
 cumprob[length(cumprob)] <- 1
 q <- ifelse(p > 0, cut(p, c(0, cumprob)) + x[1] - 1,  -Inf)

The fact that the last element of cumprob is set to 1 causes the call to
cut() to produce the error message (but I don't know exactly why). If the
expression "cumprob[ length( cumprob)] <- 1" is skipped, everything runs
smoothly (and yields the result q = 1460).

So, is that problematic expression needed at all? (Per definition cumprob[
length( cumprob)] contains 1 anyway, except if rounding errors are an
issue here ... Are they?)

 Regards  --  Gerrit

PS: S-PLUS Version 6.0 Release 1 for Sun SPARC, SunOS 5.6 :  2000.

 -------------------------------------------------------------------------
  Dr. Gerrit Eichner                        Mathematical Institute of the
  gerrit.eichner@math.uni-giessen.de        Justus-Liebig-Univ. Giessen
  Tel: +49-(0)641-99-32104            Arndtstr. 2, 35392 Giessen, Germany
  Fax: +49-(0)641-99-32029      http://www.math.uni-giessen.de/Stochastik
 -------------------------------------------------------------------------


<Prev in Thread] Current Thread [Next in Thread>
  • "Bug in qwilcox()?" revisited, Gerrit Eichner <=