================================================= This is an automated collection of open reports in the R-bugs repository. ================================================= Directory: Accuracy * PR# 1228 * Subject: bug with var(rep(1e30, 3)) From: Emmanuel Paradis Date: Wed, 26 Dec 2001 13:03:31 +0100 --needs a bit of work on src/main/cov.c ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 1228 <== From postmaster@franz.stat.wisc.edu Wed Dec 26 13:07:55 2001 Received: from franz.stat.wisc.edu (root@franz.stat.wisc.edu [128.105.174.95]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id NAA13232 for ; Wed, 26 Dec 2001 13:07:55 +0100 (MET) Received: from isem.isem.univ-montp2.fr (really [162.38.183.200]) by franz.stat.wisc.edu via smail with esmtp id (Debian Smail3.2.0.114) for ; Wed, 26 Dec 2001 06:07:27 -0600 (CST) Received: from diversi (diversi.isem.univ-montp2.fr [162.38.183.45]) by isem.isem.univ-montp2.fr (8.9.1b+Sun/8.9.1) with SMTP id MAA19343 for ; Wed, 26 Dec 2001 12:58:33 +0100 (MET) Message-Id: <3.0.32.20011226130330.00910b80@162.38.183.200> X-Sender: paradis@162.38.183.200 X-Mailer: Windows Eudora Pro Version 3.0 (32) Date: Wed, 26 Dec 2001 13:03:31 +0100 To: r-bugs@r-project.org From: Emmanuel Paradis Subject: bug with var(rep(1e30, 3)) Mime-Version: 1.0 Content-Type: text/plain; charset="iso-8859-1" Content-Transfer-Encoding: 8bit X-MIME-Autoconverted: from quoted-printable to 8bit by pubhealth.ku.dk id NAA13232 There seems to be a bug with var() when the argument is a vector with exactly three values of 1e30 (or close to this value). This does not happen with twice, four (or more) times this value, or another value. > var(rep(1e30, 3)) [1] 2.971056e+28 > var(rep(1.2e30, 3)) [1] 2.971056e+28 > var(rep(0.9e30, 3)) [1] 2.971056e+28 > var(rep(0.8e30, 3)) [1] 0 > var(rep(1e29, 3)) [1] 0 > var(rep(1e30, 2)) [1] 0 > var(rep(1e30, 4)) [1] 0 > var(rep(1e31, 3)) [1] 0 The bug is repeatable, and I got the same results with R 1.3.1 (binary from CRAN) and R 1.4.0 (binary from BD Ripley's site) both on Windows NT, and R 1.4.0 on Solaris (compiled from sources). Emmanuel Paradis Laboratoire de Paléontologie Institut des Sciences de l'Évolution Université Montpellier II F-34095 Montpellier cédex 05 France phone: +33 4 67 14 39 64 fax: +33 4 67 14 36 10 mailto:paradis@isem.univ-montp2.fr ==> 1228.reply.1 <== From p.dalgaard@biostat.ku.dk Wed Dec 26 15:32:38 2001 Received: from blueberry.kubism.ku.dk (blueberry [130.225.108.193]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id PAA13516; Wed, 26 Dec 2001 15:32:38 +0100 (MET) Received: (from pd@localhost) by blueberry.kubism.ku.dk (8.11.2/8.11.2) id fBQEVRP24017; Wed, 26 Dec 2001 15:31:27 +0100 X-Authentication-Warning: blueberry.kubism.ku.dk: pd set sender to p.dalgaard@biostat.ku.dk using -f Sender: pd@pubhealth.ku.dk To: paradis@isem.univ-montp2.fr Cc: r-devel@stat.math.ethz.ch, R-bugs@biostat.ku.dk Subject: Re: [Rd] bug with var(rep(1e30, 3)) (PR#1228) References: <200112261207.NAA13242@pubhealth.ku.dk> From: Peter Dalgaard BSA Date: 26 Dec 2001 15:31:27 +0100 In-Reply-To: <200112261207.NAA13242@pubhealth.ku.dk> Message-ID: Lines: 56 User-Agent: Gnus/5.0808 (Gnus v5.8.8) Emacs/20.7 MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii paradis@isem.univ-montp2.fr writes: > There seems to be a bug with var() when the argument is a vector with > exactly three values of 1e30 (or close to this value). This does not happen > with twice, four (or more) times this value, or another value. > > > > var(rep(1e30, 3)) > [1] 2.971056e+28 > > var(rep(1.2e30, 3)) > [1] 2.971056e+28 > > var(rep(0.9e30, 3)) > [1] 2.971056e+28 > > var(rep(0.8e30, 3)) > [1] 0 > > var(rep(1e29, 3)) > [1] 0 > > var(rep(1e30, 2)) > [1] 0 > > var(rep(1e30, 4)) > [1] 0 > > var(rep(1e31, 3)) > [1] 0 > > > The bug is repeatable, and I got the same results with R 1.3.1 (binary from > CRAN) and R 1.4.0 (binary from BD Ripley's site) both on Windows NT, and R > 1.4.0 on Solaris (compiled from sources). Not a bug, roundoff: > (1e30+1e30+1e30)/3-1e30 [1] 1.407375e+14 > 3*((1e30+1e30+1e30)/3-1e30)^2/2 [1] 2.971056e+28 Relative errors of the order 1e-16 are completely expectable in floating point arithmetic. For illustration, look at it in octal, after removing a bunch of trailing zeroes: > x<-1e30/(8^16) > structure(x, class="octmode") [1] "144762623464043165" > structure(x+x+x, class="octmode") [1] "456730272634151540" If you look carefully, you will see that a bit is lost in the process since 5+5+5 is 17 octal. -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk) FAX: (+45) 35327907 ==> 1228.followup.1 <== From ripley@stats.ox.ac.uk Wed Dec 26 16:52:52 2001 Received: from toucan.stats.ox.ac.uk (toucan.stats.ox.ac.uk [163.1.20.20]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id QAA13677; Wed, 26 Dec 2001 16:52:51 +0100 (MET) Received: from gannet.stats (IDENT:root@gannet.stats [163.1.20.127]) by toucan.stats.ox.ac.uk (8.9.0/8.9.0) with ESMTP id PAA29226; Wed, 26 Dec 2001 15:52:25 GMT Received: from localhost (ripley@localhost) by gannet.stats (8.11.3/8.11.3) with ESMTP id fBQFqPg12750; Wed, 26 Dec 2001 15:52:25 GMT X-Authentication-Warning: gannet.stats: ripley owned process doing -bs Date: Wed, 26 Dec 2001 15:52:25 +0000 (GMT) From: Prof Brian Ripley X-X-Sender: To: Peter Dalgaard BSA cc: , Subject: Re: [Rd] bug with var(rep(1e30, 3)) (PR#1228) In-Reply-To: Message-ID: MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII But there are better algorithms that are not subject to this sort of rounding error. I reckon that computing a variance this way is not the quality of algorithm one should expect from R. On 26 Dec 2001, Peter Dalgaard BSA wrote: > paradis@isem.univ-montp2.fr writes: > > > There seems to be a bug with var() when the argument is a vector with > > exactly three values of 1e30 (or close to this value). This does not happen > > with twice, four (or more) times this value, or another value. > > > > > > > var(rep(1e30, 3)) > > [1] 2.971056e+28 > > > var(rep(1.2e30, 3)) > > [1] 2.971056e+28 > > > var(rep(0.9e30, 3)) > > [1] 2.971056e+28 > > > var(rep(0.8e30, 3)) > > [1] 0 > > > var(rep(1e29, 3)) > > [1] 0 > > > var(rep(1e30, 2)) > > [1] 0 > > > var(rep(1e30, 4)) > > [1] 0 > > > var(rep(1e31, 3)) > > [1] 0 > > > > > > The bug is repeatable, and I got the same results with R 1.3.1 (binary from > > CRAN) and R 1.4.0 (binary from BD Ripley's site) both on Windows NT, and R > > 1.4.0 on Solaris (compiled from sources). > > Not a bug, roundoff: > > > (1e30+1e30+1e30)/3-1e30 > [1] 1.407375e+14 > > 3*((1e30+1e30+1e30)/3-1e30)^2/2 > [1] 2.971056e+28 > > Relative errors of the order 1e-16 are completely expectable in > floating point arithmetic. > > For illustration, look at it in octal, after removing a bunch of > trailing zeroes: > > > x<-1e30/(8^16) > > structure(x, class="octmode") > [1] "144762623464043165" > > structure(x+x+x, class="octmode") > [1] "456730272634151540" > > If you look carefully, you will see that a bit is lost in the process > since 5+5+5 is 17 octal. > > -- > O__ ---- Peter Dalgaard Blegdamsvej 3 > c/ /'_ --- Dept. of Biostatistics 2200 Cph. N > (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 > ~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk) FAX: (+45) 35327907 > -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- > r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html > Send "info", "help", or "[un]subscribe" > (in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch > _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ > -- Brian D. Ripley, ripley@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 272860 (secr) Oxford OX1 3TG, UK Fax: +44 1865 272595 ==> 1228.reply.2 <== From p.dalgaard@biostat.ku.dk Wed Dec 26 18:01:10 2001 Received: from blueberry.kubism.ku.dk (blueberry [130.225.108.193]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id SAA13867; Wed, 26 Dec 2001 18:01:10 +0100 (MET) Received: (from pd@localhost) by blueberry.kubism.ku.dk (8.11.2/8.11.2) id fBQGxxE24151; Wed, 26 Dec 2001 17:59:59 +0100 X-Authentication-Warning: blueberry.kubism.ku.dk: pd set sender to p.dalgaard@biostat.ku.dk using -f Sender: pd@pubhealth.ku.dk To: Prof Brian Ripley Cc: Peter Dalgaard BSA , , Subject: Re: [Rd] bug with var(rep(1e30, 3)) (PR#1228) References: From: Peter Dalgaard BSA Date: 26 Dec 2001 17:59:59 +0100 In-Reply-To: Message-ID: Lines: 33 User-Agent: Gnus/5.0808 (Gnus v5.8.8) Emacs/20.7 MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Prof Brian Ripley writes: > But there are better algorithms that are not subject to this sort > of rounding error. > > I reckon that computing a variance this way is not the quality of > algorithm one should expect from R. Hmm. Which algorithms are there? Surely almost any kind of provisional means subtraction could avoid the the extreme case where the mean of a bunch of identical numbers is different from all of them, but will it also improve precision in the general case? Essentially our current algorithm in cov.c is xm<- sum(x)/n v <- 1/(n-1)*sum((x-xm)^2) > > > x<-1e30/(8^16) > > > structure(x, class="octmode") > > [1] "144762623464043165" > > > structure(x+x+x, class="octmode") > > [1] "456730272634151540" > > > > If you look carefully, you will see that a bit is lost in the process > > since 5+5+5 is 17 octal. -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk) FAX: (+45) 35327907 ==> 1228.reply.3 <== From maechler@stat.math.ethz.ch Tue Apr 6 12:10:15 2004 Return-Path: X-Original-To: R-bugs@biostat.ku.dk Delivered-To: R-bugs@biostat.ku.dk Received: from localhost (slam.kubism.ku.dk [192.38.18.22]) by slim.kubism.ku.dk (Postfix) with ESMTP id 3A2F8EC53 for ; Tue, 6 Apr 2004 12:10:15 +0200 (CEST) Received: from slam.kubism.ku.dk ([127.0.0.1]) by localhost (slam [127.0.0.1]) (amavisd-new, port 10024) with ESMTP id 27402-06 for ; Tue, 6 Apr 2004 12:10:14 +0200 (CEST) Received: from hypatia.math.ethz.ch (hypatia.ethz.ch [129.132.58.23]) by slam.kubism.ku.dk (Postfix) with ESMTP for ; Tue, 6 Apr 2004 12:10:13 +0200 (CEST) Received: from lynne.ethz.ch (lynne [129.132.58.30]) by hypatia.math.ethz.ch (8.12.11/8.12.11) with ESMTP id i36AADO2017399 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Tue, 6 Apr 2004 12:10:13 +0200 Received: (from maechler@localhost) by lynne.ethz.ch (8.12.10/8.12.8/Submit) id i36AACS5026625; Tue, 6 Apr 2004 12:10:12 +0200 From: Martin Maechler MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Transfer-Encoding: 7bit Message-ID: <16498.33284.713813.746088@gargle.gargle.HOWL> Date: Tue, 6 Apr 2004 12:10:12 +0200 To: daheiser@gvn.net Cc: R-bugs@biostat.ku.dk Subject: Re: [Rd] Accuracy Bug (PR#1228), (PR#6743) In-Reply-To: <20040406022435.4A6C0EFB4@slim.kubism.ku.dk> References: <20040406022435.4A6C0EFB4@slim.kubism.ku.dk> X-Mailer: VM 7.17 under Emacs 21.3.2 Reply-To: Martin Maechler X-Virus-Scanned: by amavisd-new at pubhealth.ku.dk X-Spam-Status: No, hits=-5.9 required=5.0 tests=BAYES_01,IN_REP_TO,REFERENCES,USER_AGENT_VM version=2.55 X-Spam-Level: X-Spam-Checker-Version: SpamAssassin 2.55 (1.174.2.19-2003-05-19-exp) >>>>> "daheiser" == daheiser >>>>> on Tue, 6 Apr 2004 04:24:35 +0200 (CEST) writes: daheiser> It is an error in the algorithm. "it" being the behavior reported in bug report PR#1228 --- [too bad you didn't use the whole string "PR#1228" in your subject; if you had, no new report would have been created, and things would have gone to the proper place in the R-bugs repository ...] namely the fact that var(rep(1e30, n)) does not always (for all n) give 0. With the following function tst.var <- function(x,nset= 2:100) { for(n in nset) { v <- var(rep(x,n)) if(v != 0) cat(sprintf("%4d: %20.12g\n", n,v)) } } Actually, on my current desktop (AMD Athlon,Linux, R 1.9.0beta) the computations seem to be more often exact than they used to: Even tst.var(1e30, 2:5000) doesn't produce any output nor do a few other 'x' arguments I tried --- but the fundamental criticism is still correct, and even on the Athlon, it's easy to find cases where tst.var() *does* produce output. daheiser> It is an excellent example of how errors and daheiser> faults occur when the programmer follows the daheiser> mathematical formula exactly. Welford's algorithm daheiser> does not produce this error. It gives correct daheiser> standard deviation and variance values. Actually, after reading @article{ChaTL97, author = {Tony F. Chan and John Gregg Lewis}, title = {Computing standard deviations: accuracy}, journal = {Commun. ACM}, volume = 22, number = 9, year = 1979, issn = {0001-0782}, pages = {526--531}, doi = {http://doi.acm.org/10.1145/359146.359152}, publisher = {ACM Press}, } @article{WesD97, author = {D. H. D. West}, title = {Updating mean and variance estimates: an improved method}, journal = {Commun. ACM}, volume = 22, number = 9, year = 1979, issn = {0001-0782}, pages = {532--535}, doi = {http://doi.acm.org/10.1145/359146.359153}, publisher = {ACM Press}, } I'd conclude (from page 531) that Welford's algorithm is a bit less accurate than the (very similar) "West" version, and we (the R developers) should rather implement the latter. Maybe for R 1.9.1 -- or even later {there are even more important numerical accuracy problems I know about!} Martin Maechler http://stat.ethz.ch/~maechler/ Seminar fuer Statistik, ETH-Zentrum LEO C16 Leonhardstr. 27 ETH (Federal Inst. Technology) 8092 Zurich SWITZERLAND phone: x-41-1-632-3408 fax: ...-1228 <>< ==> 1228.audit <== Fri Dec 28 20:08:08 2001 thomas moved from incoming to Accuracy Tue Apr 6 17:46:35 2004 maechler changed notes ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 2894 * Subject: qbeta hang From: terra@gnome.org Date: Thu, 1 May 2003 22:43:59 +0200 (MET DST) --R 1.7.1 incorporates the two suggestions of the 2nd followup: -- r = sqrt(-2 * log (a)) and ...(log1p(-a) + log (qq))... -- ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 2894 <== From terra@gnome.org Thu May 1 22:44:00 2003 Received: from blueberry (blueberry [192.38.19.4]) by pubhealth.ku.dk (8.12.9/8.12.9) with ESMTP id h41Khx2b017785 for ; Thu, 1 May 2003 22:43:59 +0200 (MET DST) Date: Thu, 1 May 2003 22:43:59 +0200 (MET DST) Message-Id: <200305012043.h41Khx2b017785@pubhealth.ku.dk> From: terra@gnome.org To: R-bugs@biostat.ku.dk Subject: qbeta hang Full_Name: Morten Welinder Version: 1.6.1 OS: Solaris/sparc Submission from: (NULL) (65.213.85.144) qbeta(0.1, 1e-8, 0.5, TRUE, FALSE) seems to hang for me. ==> 2894.reply.1 <== From p.dalgaard@biostat.ku.dk Fri May 2 00:16:13 2003 Received: from blueberry.kubism.ku.dk (blueberry [192.38.19.4]) by pubhealth.ku.dk (8.12.9/8.12.9) with ESMTP id h41MGD2b018113; Fri, 2 May 2003 00:16:13 +0200 (MET DST) Received: from blueberry.kubism.ku.dk (localhost.localdomain [127.0.0.1]) by blueberry.kubism.ku.dk (8.12.8/8.12.8) with ESMTP id h41MK78e014746; Fri, 2 May 2003 00:20:07 +0200 Received: (from pd@localhost) by blueberry.kubism.ku.dk (8.12.8/8.12.8/Submit) id h41MK6Ja014742; Fri, 2 May 2003 00:20:06 +0200 X-Authentication-Warning: blueberry.kubism.ku.dk: pd set sender to p.dalgaard@biostat.ku.dk using -f To: terra@gnome.org Cc: r-devel@stat.math.ethz.ch, R-bugs@biostat.ku.dk Subject: Re: [Rd] qbeta hang (PR#2894) References: <200305012044.h41Ki02b017789@pubhealth.ku.dk> From: Peter Dalgaard BSA Date: 02 May 2003 00:20:06 +0200 In-Reply-To: <200305012044.h41Ki02b017789@pubhealth.ku.dk> Message-ID: Lines: 45 User-Agent: Gnus/5.09 (Gnus v5.9.0) Emacs/21.2 MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii terra@gnome.org writes: > Full_Name: Morten Welinder > Version: 1.6.1 > OS: Solaris/sparc > Submission from: (NULL) (65.213.85.144) > > > qbeta(0.1, 1e-8, 0.5, TRUE, FALSE) seems to hang for me. confirmed on 1.7.0 Solaris 9, gcc 3.0.3 (standard build, so -O2, I assume) Morten: the gcc version is often crucial in these cases. However, the exact same thing is happening on Linux. The immediate cause is that n = fmax2(lneps/log(y), 4.0) gets large when y is in the vicinity of 1-1e-8, so the loop in src/nmath/pbeta.c:101 gets a rather high count. The algorithm isn't really stuck, it just takes a very long time. On the fastest machine that I have available: > system.time(qbeta(0.1, 1e-8, 0.5, TRUE, FALSE)) [1] 75.58 0.00 75.58 0.00 0.00 It's not really that surprising: > pbeta(1e-5, 1e-8, 0.5,, TRUE, FALSE) [1] 0.9999999 > pbeta(1e-10, 1e-8, 0.5,, TRUE, FALSE) [1] 0.9999998 > pbeta(1e-200, 1e-8, 0.5,, TRUE, FALSE) [1] 0.9999954 > pbeta(1e-309, 1e-8, 0.5,, TRUE, FALSE) [1] 0.999993 > pbeta(1e-400, 1e-8, 0.5,, TRUE, FALSE) [1] 0 and you're trying to solve pbeta(x, 1e-8, 0.5,, TRUE, FALSE) == 0.1 -- O__ ---- Peter Dalgaard Blegdamsvej 3 c/ /'_ --- Dept. of Biostatistics 2200 Cph. N (*) \(*) -- University of Copenhagen Denmark Ph: (+45) 35327918 ~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk) FAX: (+45) 35327907 ==> 2894.followup.1 <== From dmurdoch@pair.com Fri May 2 00:38:53 2003 Received: from tomts13-srv.bellnexxia.net (tomts13.bellnexxia.net [209.226.175.34]) by pubhealth.ku.dk (8.12.9/8.12.9) with ESMTP id h41Mcp2b018208 for ; Fri, 2 May 2003 00:38:52 +0200 (MET DST) Received: from djmfuj ([65.93.70.143]) by tomts13-srv.bellnexxia.net (InterMail vM.5.01.04.19 201-253-122-122-119-20020516) with SMTP id <20030501223850.JWNE22698.tomts13-srv.bellnexxia.net@djmfuj>; Thu, 1 May 2003 18:38:50 -0400 From: Duncan Murdoch To: terra@gnome.org Cc: R-bugs@biostat.ku.dk Subject: Re: [Rd] qbeta hang (PR#2894) Date: Thu, 01 May 2003 18:38:49 -0400 Message-ID: References: <200305012044.h41Ki02b017789@pubhealth.ku.dk> In-Reply-To: <200305012044.h41Ki02b017789@pubhealth.ku.dk> X-Mailer: Forte Agent 1.91/32.564 MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Transfer-Encoding: 7bit On Thu, 1 May 2003 22:44:00 +0200 (MET DST), you wrote: >Full_Name: Morten Welinder >Version: 1.6.1 >OS: Solaris/sparc >Submission from: (NULL) (65.213.85.144) > > >qbeta(0.1, 1e-8, 0.5, TRUE, FALSE) seems to hang for me. In 1.7.0-patched in Windows, it's very, very slow, but eventually produces an answer: > qbeta(0.1, 1e-8, 0.5, TRUE, FALSE) [1] 4.400784e-309 The answer should be really small (the mean is 2e-8, and it's strongly skewed to the right), so I'd guess it's giving the right answer, and is just slow to converge to it. Duncan Murdoch ==> 2894.followup.2 <== From welinder@rentec.com Fri May 2 16:03:42 2003 Received: from ram.rentec.com (ram.rentec.com [192.5.35.66]) by pubhealth.ku.dk (8.12.9/8.12.9) with ESMTP id h42E3e2b026336; Fri, 2 May 2003 16:03:41 +0200 (MET DST) Received: from robbie.rentec.com (robbie.rentec.com [192.5.35.99]) by ram.rentec.com (8.12.1/8.12.1) with ESMTP id h42E3OR5024597; Fri, 2 May 2003 10:03:24 -0400 (EDT) Received: (from welinder@localhost) by robbie.rentec.com (8.12.5+Sun/8.12.5/Submit) id h42E3NcU001790; Fri, 2 May 2003 10:03:23 -0400 (EDT) Date: Fri, 2 May 2003 10:03:23 -0400 (EDT) Message-Id: <200305021403.h42E3NcU001790@robbie.rentec.com> From: Morten Welinder To: p.dalgaard@biostat.ku.dk CC: r-devel@stat.math.ethz.ch, R-bugs@biostat.ku.dk In-reply-to: (message from Peter Dalgaard BSA on 02 May 2003 00:20:06 +0200) Subject: Re: [Rd] qbeta hang (PR#2894) References: <200305012044.h41Ki02b017789@pubhealth.ku.dk> Ok, I can confirm that it does not, in fact, loop forever. Just a close approximation. > Morten: the gcc version is often crucial in these cases. Sorry. gcc-2.7.2 (both -O2 and -O0) on Solaris 2.9 / sparc. The problems seem to be... 1. The initial guess underflows to zero. 2. That guess is replaced by the truly awful guess 0.5. 3. The root finding algorithm ignores all the nice properties of pbeta -- such as it being monotonically increasing and bounded. > and you're trying to solve pbeta(x, 1e-8, 0.5,, TRUE, FALSE) == 0.1 Not quite. What I was trying to do was the see what would happen with an insanely small alpha. Actually doing that would involve getting the arguments right, of course, :-) The not-quite-hang was an unplesant side effect of swapping the args. I notice things like r = sqrt(-log(a * a)); in the code. I fail to see why that isn't just r = sqrt(-2 * log (a)) which ought to be faster and more accurate. Then later ...log((1. - a) * qq)... which might be better as ...(log1p(-a) + log (qq))... if a is close to zero. There are lots of other places that worry me with respect to cancellation errors, for example r = 1 - pp; t = 1 - qq; Morten ==> 2894.audit <== Thu May 22 23:01:16 2003 thomas moved from incoming to Accuracy Mon Jun 16 11:24:33 2003 maechler changed notes Mon Jun 16 11:24:33 2003 maechler foobar Mon Jun 16 11:24:56 2003 maechler changed notes Mon Jun 16 11:24:56 2003 maechler foobar ==> 2894.followup.3 <== From tim.massingham@ebi.ac.uk Mon Sep 22 15:47:29 2003 Received: from maui.ebi.ac.uk (maui.ebi.ac.uk [193.62.196.100]) by pubhealth.ku.dk (8.12.9/8.12.9) with ESMTP id h8MDlSJH021438 for ; Mon, 22 Sep 2003 15:47:29 +0200 (MET DST) Received: from teapot.ebi.ac.uk (teapot.ebi.ac.uk [172.22.100.116]) by maui.ebi.ac.uk (8.11.7+Sun/8.11.7) with ESMTP id h8MDlJm01253 for ; Mon, 22 Sep 2003 14:47:19 +0100 (BST) From: Tim Massingham Organization: EBI Subject: PR#2894 Date: Mon, 22 Sep 2003 14:47:54 +0100 User-Agent: KMail/1.5.1 To: r-bugs@biostat.ku.dk MIME-Version: 1.0 Content-Type: text/plain; charset="iso-8859-1" Content-Transfer-Encoding: 7bit Content-Disposition: inline Message-Id: <200309221447.54338.tim.massingham@ebi.ac.uk> X-EBI-Information: This email is scanned using www.mailscanner.info. X-EBI: Found to be clean X-EBI-SpamCheck: not spam, SpamAssassin (score=-5.8, required 5, USER_AGENT_KMAIL -5.80) >Date: Fri, 2 May 2003 10:03:23 -0400 (EDT) From: Morten Welinder >To: p.dalgaard@biostat.ku.dk >CC: r-devel@stat.math.ethz.ch, R-bugs@biostat.ku.dk >Subject: Re: [Rd] qbeta hang (PR#2894) > >Ok, I can confirm that it does not, in fact, loop forever. Just a close >approximation. ... >There are lots of other places that worry me with respect to cancellation >errors, for example > > r = 1 - pp; > t = 1 - qq; These can also be removed by changing qbeta.c:156-157 (R-1.7.1) to y = (y-a) * xinbta * (1.0-xinbta) * exp (logbeta - pp * log(xinbta) - qq * log1p(-xinbta)); I don't understand why you have qbeta.c:167 if (fabs(y)<=acu) goto L_converged which will fail for certain values of p & q, when x is close to zero as the beta density tends to infinity. Why not use an exit condition based on how far away from 'a' you are? qbeta.c:174 (prevent excess precision) Might this not just get optimized out? I doubt modern compilers respect the volatile keyword. You should probably replace this with a safe comparison. if(fabs(tx-xinbta)<=DBL_EPSILON*fabs(xinbta))... ** This should be checked with someone familiar with floating point arithmetic; I'm not and may have got the correct expression wrong. The 'infinite loop' is due to the speed of the pbeta routine, rather than the initial approximation from the qbeta. There is another algorithm available for pbeta (Algorithm 708, CACM 18. 360--373, 1992) which may be of use here. This continued fraction approximation is recommended by Numerical Recipes over the power series approximation (make of that what you will!) TimM. ==> 2894.followup.4 <== From tim.massingham@ebi.ac.uk Mon Sep 22 15:49:10 2003 Received: from maui.ebi.ac.uk (maui.ebi.ac.uk [193.62.196.100]) by pubhealth.ku.dk (8.12.9/8.12.9) with ESMTP id h8MDn9JH021470 for ; Mon, 22 Sep 2003 15:49:10 +0200 (MET DST) Received: from teapot.ebi.ac.uk (teapot.ebi.ac.uk [172.22.100.116]) by maui.ebi.ac.uk (8.11.7+Sun/8.11.7) with ESMTP id h8MDn4m01806 for ; Mon, 22 Sep 2003 14:49:04 +0100 (BST) From: Tim Massingham Organization: EBI Subject: Re: PR#2894 Date: Mon, 22 Sep 2003 14:49:40 +0100 User-Agent: KMail/1.5.1 To: r-bugs@biostat.ku.dk MIME-Version: 1.0 Content-Type: text/plain; charset="iso-8859-1" Content-Transfer-Encoding: 7bit Content-Disposition: inline Message-Id: <200309221449.40171.tim.massingham@ebi.ac.uk> X-EBI-Information: This email is scanned using www.mailscanner.info. X-EBI: Found to be clean X-EBI-SpamCheck: not spam, SpamAssassin (score=-15.5, required 5, EMAIL_ATTRIBUTION -6.50, QUOTED_EMAIL_TEXT -3.20, USER_AGENT_KMAIL -5.80) On Monday 22 Sep 2003 2:36 pm, I wrote: > The 'infinite loop' is due to the speed of the pbeta routine, rather > than the initial approximation from the qbeta. There is another > algorithm available for pbeta (Algorithm 708, CACM 18. 360--373, 1992) That should be TOMS, not CACM. Sorry. > which may be of use here. This continued fraction approximation is > recommended by Numerical Recipes over the power series approximation > (make of that what you will!) TimM. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 5727 * Subject: Bug in pbinom? From: Volkmar Liebscher Date: Fri, 12 Dec 2003 11:18:12 +0100 --(really on of the infamous border case problem of pbeta) ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 5727 <== From mailnull@stat.math.ethz.ch Fri Dec 12 11:19:22 2003 Return-Path: X-Original-To: R-bugs@biostat.ku.dk Delivered-To: R-bugs@biostat.ku.dk Received: from hypatia.math.ethz.ch (hypatia.ethz.ch [129.132.58.23]) by slim.kubism.ku.dk (Postfix) with ESMTP id 6EDB2FC00 for ; Fri, 12 Dec 2003 11:19:21 +0100 (CET) Received: from hypatia.math.ethz.ch (localhost [127.0.0.1]) by hypatia.math.ethz.ch (8.12.10/8.12.10) with ESMTP id hBCAHamw004986 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Fri, 12 Dec 2003 11:17:36 +0100 (MET) Received: (from mailnull@localhost) by hypatia.math.ethz.ch (8.12.10/8.12.10/Submit) id hBCAHZ7F004985 for R-bugs@biostat.ku.dk; Fri, 12 Dec 2003 11:17:35 +0100 (MET) Received: from franz.stat.wisc.edu (www.omegahat.org [128.105.174.32]) by hypatia.math.ethz.ch (8.12.10/8.12.10) with ESMTP id hBCAHNmv004962 for ; Fri, 12 Dec 2003 11:17:25 +0100 (MET) Received: from servus.gsf.de ([146.107.101.11]) by franz.stat.wisc.edu with esmtp (Exim 3.35 #1 (Debian)) id 1AUkO0-0001lI-00 for ; Fri, 12 Dec 2003 04:19:04 -0600 Received: from hausdorff.gsf.de (hausdorff.gsf.de [146.107.58.25]) by servus.gsf.de (8.8.6 (PHNE_17135)/8.8.6.+++) with ESMTP id LAA09567 for ; Fri, 12 Dec 2003 11:21:36 +0100 (MET) From: Volkmar Liebscher Reply-To: liebscher@gsf.de Organization: GSF To: r-bugs@r-project.org Subject: Bug in pbinom? Date: Fri, 12 Dec 2003 11:18:12 +0100 User-Agent: KMail/1.5.1 MIME-Version: 1.0 Content-Type: text/plain; charset="iso-8859-1" Content-Disposition: inline Message-Id: <200312121118.12004.liebscher@gsf.de> X-MIME-Autoconverted: from 8bit to quoted-printable by servus.gsf.de id LAA09567 X-Virus-Scanned: by amavisd-milter (http://amavis.org/) X-Virus-Scanned: by amavisd-milter (http://amavis.org/) Content-Transfer-Encoding: 8bit X-MIME-Autoconverted: from quoted-printable to 8bit by hypatia.math.ethz.ch id hBCAHNmv004962 X-Spam-Checker-Version: SpamAssassin 2.61 (1.212.2.1-2003-12-09-exp) on hypatia.math.ethz.ch X-Spam-Status: No, hits=-4.9 required=5.0 tests=BAYES_00 autolearn=ham version=2.61 X-Spam-Level: Dear colleagues, the command pbinom(size=1.95*10^10,prob=seq(1.02,1.03,10^-3)*10^-10,q=1,lower=F,log=T) produced the output [1] -46.0120973 -46.0101369 -46.0081784 -46.0062239 -46.0042692 -46.0023165 [7] -0.5205666 -0.5196801 -0.5187944 -0.5179102 -0.5170278 what seems strange. Has this survived in newer versions? Best regards, Volkmar Liebscher -- Dr. Volkmar Liebscher Institute of Biomathematics and Biometry GSF - National Research Centre for Environment and Health Ingolstädter Landstr.1 D--85758 Neuherberg Germany E-mail: liebscher@gsf.de URL: http://www.gsf.de/ibb/homepages/liebscher phone: +49 89 3187 2907 Fax: +49 89 3187 3369 --please do not edit the information below-- Version: platform = i686-pc-linux-gnu arch = i686 os = linux-gnu system = i686, linux-gnu status = major = 1 minor = 7.1 year = 2003 month = 06 day = 16 language = R Search Path: .GlobalEnv, package:methods, package:ctest, package:mva, package:modreg, package:nls, package:ts, Autoloads, package:base ==> 5727.followup.1 <== From ripley@stats.ox.ac.uk Fri Dec 12 11:44:59 2003 Return-Path: X-Original-To: R-bugs@biostat.ku.dk Delivered-To: R-bugs@biostat.ku.dk Received: from toucan.stats.ox.ac.uk (toucan.stats.ox.ac.uk [163.1.20.20]) by slim.kubism.ku.dk (Postfix) with ESMTP id 0991FFC00 for ; Fri, 12 Dec 2003 11:44:59 +0100 (CET) Received: from gannet.stats (gannet.stats [163.1.20.127]) by toucan.stats.ox.ac.uk (8.12.10/8.12.10) with ESMTP id hBCAivoO002696; Fri, 12 Dec 2003 10:44:57 GMT Date: Fri, 12 Dec 2003 10:44:43 +0000 (GMT) From: Prof Brian Ripley X-X-Sender: ripley@gannet.stats To: liebscher@gsf.de Cc: r-devel@stat.math.ethz.ch, Subject: Re: [Rd] Bug in pbinom? (PR#5727) In-Reply-To: <20031212101923.309D2FC01@slim.kubism.ku.dk> Message-ID: MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII Oh, come on: have you never heard of the Poisson approximation? > ppois(1, 1.95*seq(1.02,1.03,10^-3), lower=F,log=T) [1] -0.5259247 -0.5250276 -0.5241321 -0.5232383 -0.5223462 -0.5214557 [7] -0.5205669 -0.5196798 -0.5187943 -0.5179104 -0.5170282 Probably you believe that computers never make errors too. On Fri, 12 Dec 2003 liebscher@gsf.de wrote: > Dear colleagues, > > the command > > pbinom(size=1.95*10^10,prob=seq(1.02,1.03,10^-3)*10^-10,q=1,lower=F,log=T) > > produced the output > > [1] -46.0120973 -46.0101369 -46.0081784 -46.0062239 -46.0042692 -46.0023165 > [7] -0.5205666 -0.5196801 -0.5187944 -0.5179102 -0.5170278 > > > what seems strange. Has this survived in newer versions? > > Best regards, > Volkmar Liebscher > -- Brian D. Ripley, ripley@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 ==> 5727.audit <== Mon Dec 15 08:29:22 2003 ripley moved from incoming to Accuracy Tue Apr 6 17:44:55 2004 maechler changed notes ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 6196 * Subject: Accuracy: Correct sums in rowSums(), colSums() From: Nick.Efthymiou@Schwab.com Date: Tue, 30 Dec 2003 01:57:19 +0100 (CET) --no requested patch against the current sources forthcoming. -- --Seems unlikely to be worthwhile for rowSums and not sum! ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 6196 <== From Nick.Efthymiou@Schwab.com Tue Dec 30 01:57:20 2003 Return-Path: X-Original-To: R-bugs@biostat.ku.dk Delivered-To: R-bugs@biostat.ku.dk Received: from blueberry (blueberry.kubism.ku.dk [192.38.19.4]) by slim.kubism.ku.dk (Postfix) with ESMTP id F26811044C for ; Tue, 30 Dec 2003 01:57:19 +0100 (CET) From: Nick.Efthymiou@Schwab.com To: R-bugs@biostat.ku.dk Subject: Accuracy: Correct sums in rowSums(), colSums() Message-Id: <20031230005719.F26811044C@slim.kubism.ku.dk> Date: Tue, 30 Dec 2003 01:57:19 +0100 (CET) Full_Name: Nick Efthymiou Version: R1.5.0 and above OS: Red Hat Linux Submission from: (NULL) (162.93.14.73) With the introduction of the functions rowSums(), colSums(), rowMeans() and colMeans() in R1.5.0, function "SEXP do_colsum(SEXP call, SEXP op, SEXP args, SEXP rho)" was added to perform the fast summations. We have an excellent opportunity to improve the accuracy by implementing Kahan summation here. Kahan summation is described in @article{ goldberg91what, author = "David Goldberg", title = "What Every Computer Scientist Should Know About Floating-Point Arithmetic", journal = "ACM Computing Surveys", volume = "23", number = "1", pages = "5--48", year = "1991", url = "citeseer.nj.nec.com/goldberg91what.html" } (http://citeseer.nj.nec.com/goldberg91what.html) and in the article "Floating-point Summation", by Evan Manning, in C/C++ User's Journal, Sept. 1996. The proposed improvement has been tested in !IEEE_754 mode, the patch is attached. It is intended to be applied against R-1.7.1/src/main/array.c --------- Cut here ---------- *** array.c.old Mon Dec 15 17:33:23 2003 --- array.c Mon Dec 15 17:33:45 2003 *************** *** 1016,1022 **** int OP, n, p, cnt = 0, i, j, type; Rboolean NaRm, keepNA; int *ix; ! double *rx, sum = 0.0; checkArity(op, args); x = CAR(args); args = CDR(args); --- 1016,1022 ---- int OP, n, p, cnt = 0, i, j, type; Rboolean NaRm, keepNA; int *ix; ! double *rx, sum = 0.0, correction = 0.0; checkArity(op, args); x = CAR(args); args = CDR(args); *************** *** 1046,1062 **** switch (type) { case REALSXP: rx = REAL(x) + n*j; #ifdef IEEE_754 if (keepNA) ! for (sum = 0., i = 0; i < n; i++) sum += *rx++; else { ! for (cnt = 0, sum = 0., i = 0; i < n; i++, rx++) ! if (!ISNAN(*rx)) {cnt++; sum += *rx;} else if (keepNA) {sum = NA_REAL; break;} } #else ! for (cnt = 0, sum = 0., i = 0; i < n; i++, rx++) ! if (!ISNAN(*rx)) {cnt++; sum += *rx;} else if (keepNA) {sum = NA_REAL; break;} #endif break; --- 1046,1082 ---- switch (type) { case REALSXP: rx = REAL(x) + n*j; + sum = 0.; #ifdef IEEE_754 if (keepNA) ! for (i = 0; i < n; i++) sum += *rx++; else { ! for (cnt = 0, i = 0; i < n; i++, rx++) ! if (!ISNAN(*rx)) { ! /* TODO: Kahan summation */ ! cnt++; sum += *rx; ! } else if (keepNA) {sum = NA_REAL; break;} } #else ! i = cnt = 0; ! if (i < n) { ! if (ISNAN(*rx)) { ! if (keepNA) {sum = NA_REAL; break /* switch */;} ! } else { ! cnt++; sum = *rx; ! } ! i++; rx++; ! } ! for (; i < n; i++, rx++) ! if (!ISNAN(*rx)) { ! /* Kahan summation */ ! double yhilo = *rx - correction; ! double tsum = sum + yhilo; ! correction = (tsum - sum) - yhilo; ! sum = tsum; ! cnt++; ! } else if (keepNA) {sum = NA_REAL; break;} #endif break; *************** *** 1121,1137 **** switch (type) { case REALSXP: rx = REAL(x) + i; #ifdef IEEE_754 if (keepNA) ! for (sum = 0., j = 0; j < p; j++, rx += n) sum += *rx; else { ! for (cnt = 0, sum = 0., j = 0; j < p; j++, rx += n) ! if (!ISNAN(*rx)) {cnt++; sum += *rx;} else if (keepNA) {sum = NA_REAL; break;} } #else ! for (cnt = 0, sum = 0., j = 0; j < p; j++, rx += n) ! if (!ISNAN(*rx)) {cnt++; sum += *rx;} else if (keepNA) {sum = NA_REAL; break;} #endif break; --- 1141,1177 ---- switch (type) { case REALSXP: rx = REAL(x) + i; + sum = 0.; #ifdef IEEE_754 if (keepNA) ! for (j = 0; j < p; j++, rx += n) sum += *rx; else { ! for (cnt = 0, j = 0; j < p; j++, rx += n) ! if (!ISNAN(*rx)) { ! /* TODO: Kahan summation */ ! cnt++; sum += *rx; ! } else if (keepNA) {sum = NA_REAL; break;} } #else ! j = cnt = 0; ! if (j < p) { ! if (ISNAN(*rx)) { ! if (keepNA) {sum = NA_REAL; break /* switch */;} ! } else { ! cnt++; sum = *rx; ! } ! j++; rx += n; ! } ! for (; j < p; j++, rx += n) ! if (!ISNAN(*rx)) { ! /* Kahan summation */ ! double yhilo = *rx - correction; ! double tsum = sum + yhilo; ! correction = (tsum - sum) - yhilo; ! sum = tsum; ! cnt++; ! } else if (keepNA) {sum = NA_REAL; break;} #endif break; ---------end of patch ------- Since someone is certainly wondering why I bothered to bring this up - I was trying to optimize the order of a Chebyshev approximation to a specific function so that I accomplish the most power-efficient calculation at one ulp accuracy (multiplications and additions consume battery power). The optimization was based on minimizing the relative error between the true function and the approximation. It involved a lot of sums and was not converging properly. With the patch, I got the convergence I needed. Thanks, - Nick - ==> 6196.followup.1 <== From ripley@stats.ox.ac.uk Tue Dec 30 09:30:23 2003 Return-Path: X-Original-To: R-bugs@biostat.ku.dk Delivered-To: R-bugs@biostat.ku.dk Received: from toucan.stats.ox.ac.uk (toucan.stats.ox.ac.uk [163.1.20.20]) by slim.kubism.ku.dk (Postfix) with ESMTP id C64B710457 for ; Tue, 30 Dec 2003 09:30:22 +0100 (CET) Received: from gannet.stats (gannet.stats [163.1.20.127]) by toucan.stats.ox.ac.uk (8.12.10/8.12.10) with ESMTP id hBU8ULoO014755; Tue, 30 Dec 2003 08:30:21 GMT Date: Tue, 30 Dec 2003 08:30:21 +0000 (GMT) From: Prof Brian Ripley X-X-Sender: ripley@gannet.stats To: Nick.Efthymiou@schwab.com Cc: r-devel@stat.math.ethz.ch, Subject: Re: [Rd] Accuracy: Correct sums in rowSums(), colSums() (PR#6196) In-Reply-To: <20031230005720.9CEC71044D@slim.kubism.ku.dk> Message-ID: MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII Could you please prepare a patch against the current R-devel sources? 1.7.1 is several versions old. I was puzzled as to why you are recommending improving the accuracy of the `fast' methods and not the slow one, applying apply() to mean() and sum(). The improvement might be larger in those cases: what Goldberg calls `a dramatic improvement' is from N eps down to 2 eps, and that is for the worst case and is only dramatic if N is dramatically bigger than 2. It is also not clear that it is effective on machines with extended-precision registers unless N is very large. Some of us are well aware of the issues of floating point arithmetic, but R is normally used for statistical computations, where errors in the data (rounding, measurement etc) normally are orders of magnitude bigger than those in floating-point arithmetic. So implementation effort has been spent on handling the statistical errors first. On Tue, 30 Dec 2003 Nick.Efthymiou@schwab.com wrote: > Full_Name: Nick Efthymiou > Version: R1.5.0 and above > OS: Red Hat Linux > Submission from: (NULL) (162.93.14.73) > > > With the introduction of the functions rowSums(), colSums(), rowMeans() and > colMeans() in R1.5.0, function "SEXP do_colsum(SEXP call, SEXP op, SEXP args, > SEXP rho)" was added to perform the fast summations. We have an excellent > opportunity to improve the accuracy by implementing Kahan summation here. > > Kahan summation is described in > > @article{ goldberg91what, > author = "David Goldberg", > title = "What Every Computer Scientist Should Know About Floating-Point > Arithmetic", > journal = "ACM Computing Surveys", > volume = "23", > number = "1", > pages = "5--48", > year = "1991", > url = "citeseer.nj.nec.com/goldberg91what.html" } > > > (http://citeseer.nj.nec.com/goldberg91what.html) > > and in the article "Floating-point Summation", by Evan Manning, in C/C++ User's > Journal, Sept. 1996. > > The proposed improvement has been tested in !IEEE_754 mode, the patch is > attached. > It is intended to be applied against R-1.7.1/src/main/array.c > > --------- Cut here ---------- > *** array.c.old Mon Dec 15 17:33:23 2003 > --- array.c Mon Dec 15 17:33:45 2003 > *************** > *** 1016,1022 **** > int OP, n, p, cnt = 0, i, j, type; > Rboolean NaRm, keepNA; > int *ix; > ! double *rx, sum = 0.0; > > checkArity(op, args); > x = CAR(args); args = CDR(args); > --- 1016,1022 ---- > int OP, n, p, cnt = 0, i, j, type; > Rboolean NaRm, keepNA; > int *ix; > ! double *rx, sum = 0.0, correction = 0.0; > > checkArity(op, args); > x = CAR(args); args = CDR(args); > *************** > *** 1046,1062 **** > switch (type) { > case REALSXP: > rx = REAL(x) + n*j; > #ifdef IEEE_754 > if (keepNA) > ! for (sum = 0., i = 0; i < n; i++) sum += *rx++; > else { > ! for (cnt = 0, sum = 0., i = 0; i < n; i++, rx++) > ! if (!ISNAN(*rx)) {cnt++; sum += *rx;} > else if (keepNA) {sum = NA_REAL; break;} > } > #else > ! for (cnt = 0, sum = 0., i = 0; i < n; i++, rx++) > ! if (!ISNAN(*rx)) {cnt++; sum += *rx;} > else if (keepNA) {sum = NA_REAL; break;} > #endif > break; > --- 1046,1082 ---- > switch (type) { > case REALSXP: > rx = REAL(x) + n*j; > + sum = 0.; > #ifdef IEEE_754 > if (keepNA) > ! for (i = 0; i < n; i++) sum += *rx++; > else { > ! for (cnt = 0, i = 0; i < n; i++, rx++) > ! if (!ISNAN(*rx)) { > ! /* TODO: Kahan summation */ > ! cnt++; sum += *rx; > ! } > else if (keepNA) {sum = NA_REAL; break;} > } > #else > ! i = cnt = 0; > ! if (i < n) { > ! if (ISNAN(*rx)) { > ! if (keepNA) {sum = NA_REAL; break /* switch */;} > ! } else { > ! cnt++; sum = *rx; > ! } > ! i++; rx++; > ! } > ! for (; i < n; i++, rx++) > ! if (!ISNAN(*rx)) { > ! /* Kahan summation */ > ! double yhilo = *rx - correction; > ! double tsum = sum + yhilo; > ! correction = (tsum - sum) - yhilo; > ! sum = tsum; > ! cnt++; > ! } > else if (keepNA) {sum = NA_REAL; break;} > #endif > break; > *************** > *** 1121,1137 **** > switch (type) { > case REALSXP: > rx = REAL(x) + i; > #ifdef IEEE_754 > if (keepNA) > ! for (sum = 0., j = 0; j < p; j++, rx += n) sum += *rx; > else { > ! for (cnt = 0, sum = 0., j = 0; j < p; j++, rx += n) > ! if (!ISNAN(*rx)) {cnt++; sum += *rx;} > else if (keepNA) {sum = NA_REAL; break;} > } > #else > ! for (cnt = 0, sum = 0., j = 0; j < p; j++, rx += n) > ! if (!ISNAN(*rx)) {cnt++; sum += *rx;} > else if (keepNA) {sum = NA_REAL; break;} > #endif > break; > --- 1141,1177 ---- > switch (type) { > case REALSXP: > rx = REAL(x) + i; > + sum = 0.; > #ifdef IEEE_754 > if (keepNA) > ! for (j = 0; j < p; j++, rx += n) sum += *rx; > else { > ! for (cnt = 0, j = 0; j < p; j++, rx += n) > ! if (!ISNAN(*rx)) { > ! /* TODO: Kahan summation */ > ! cnt++; sum += *rx; > ! } > else if (keepNA) {sum = NA_REAL; break;} > } > #else > ! j = cnt = 0; > ! if (j < p) { > ! if (ISNAN(*rx)) { > ! if (keepNA) {sum = NA_REAL; break /* switch */;} > ! } else { > ! cnt++; sum = *rx; > ! } > ! j++; rx += n; > ! } > ! for (; j < p; j++, rx += n) > ! if (!ISNAN(*rx)) { > ! /* Kahan summation */ > ! double yhilo = *rx - correction; > ! double tsum = sum + yhilo; > ! correction = (tsum - sum) - yhilo; > ! sum = tsum; > ! cnt++; > ! } > else if (keepNA) {sum = NA_REAL; break;} > #endif > break; > ---------end of patch ------- > > > Since someone is certainly wondering why I bothered to bring this up - I was > trying to optimize the order of a Chebyshev approximation to a specific function > so that I accomplish the most power-efficient calculation at one ulp accuracy > (multiplications and additions consume battery power). The optimization was > based on minimizing the relative error between the true function and the > approximation. It involved a lot of sums and was not converging properly. With > the patch, I got the convergence I needed. > > Thanks, > > - Nick - > > ______________________________________________ > R-devel@stat.math.ethz.ch mailing list > https://www.stat.math.ethz.ch/mailman/listinfo/r-devel > > -- Brian D. Ripley, ripley@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 ==> 6196.audit <== Thu Jan 8 15:30:46 2004 ripley moved from incoming to Accuracy Mon Mar 15 21:33:24 2004 ripley changed notes ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 6743 * Subject: Accuracy Bug 1228 From: "David Heiser" Date: Mon, 5 Apr 2004 19:23:35 -0700 --This is really bug PR#1228 (in category 'Accuracy') ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 6743 <== From daheiser@gvn.net Tue Apr 6 04:24:06 2004 Return-Path: X-Original-To: r-bugs@biostat.ku.dk Delivered-To: r-bugs@biostat.ku.dk Received: from localhost (slam.kubism.ku.dk [192.38.18.22]) by slim.kubism.ku.dk (Postfix) with ESMTP id 9DE91EDF1 for ; Tue, 6 Apr 2004 04:24:06 +0200 (CEST) Received: from slam.kubism.ku.dk ([127.0.0.1]) by localhost (slam [127.0.0.1]) (amavisd-new, port 10024) with ESMTP id 24016-02 for ; Tue, 6 Apr 2004 04:24:05 +0200 (CEST) Received: from rook.innercite.com (rook.innercite.com [158.222.5.8]) by slam.kubism.ku.dk (Postfix) with ESMTP for ; Tue, 6 Apr 2004 04:23:38 +0200 (CEST) Received: from f0q5g8 (host-224-114.dialup.innercite.com [158.222.224.114]) by rook.innercite.com (8.12.3/8.12.3/Debian-6.4) with SMTP id i362NZ3N018444 for ; Mon, 5 Apr 2004 19:23:36 -0700 From: "David Heiser" To: Subject: Accuracy Bug 1228 Date: Mon, 5 Apr 2004 19:23:35 -0700 Message-ID: MIME-Version: 1.0 Content-Type: text/plain; charset="Windows-1252" Content-Transfer-Encoding: 7bit X-Priority: 3 (Normal) X-MSMail-Priority: Normal X-Mailer: Microsoft Outlook IMO, Build 9.0.6604 (9.0.2911.0) Importance: Normal X-MimeOLE: Produced By Microsoft MimeOLE V6.00.2800.1165 X-Virus-Scanned: by amavisd-new at pubhealth.ku.dk X-Spam-Status: No, hits=-0.1 required=5.0 tests=MSGID_GOOD_EXCHANGE version=2.55 X-Spam-Level: X-Spam-Checker-Version: SpamAssassin 2.55 (1.174.2.19-2003-05-19-exp) It is an error in the algorithm. It is an excellent example of how errors and faults occur when the programmer follows the mathematical formula exactly. Welford's algorithm does not produce this error. It gives correct standard deviation and variance values. David Heiser ==> 6743.audit <== Tue Apr 6 17:40:49 2004 maechler changed notes Tue Apr 6 17:40:49 2004 maechler moved from incoming to Accuracy ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 6772 * Subject: phyper accuracy and efficiency From: terra@gnome.org Date: Thu, 15 Apr 2004 18:06:08 +0200 (CEST) -- no notes -- ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 6772 <== From terra@gnome.org Thu Apr 15 18:06:08 2004 Return-Path: X-Original-To: R-bugs@biostat.ku.dk Delivered-To: R-bugs@biostat.ku.dk Received: from blueberry (blueberry.kubism.ku.dk [192.38.19.4]) by slim.kubism.ku.dk (Postfix) with ESMTP id 2D24CEFC2 for ; Thu, 15 Apr 2004 18:06:08 +0200 (CEST) From: terra@gnome.org To: R-bugs@biostat.ku.dk Subject: phyper accuracy and efficiency Message-Id: <20040415160608.2D24CEFC2@slim.kubism.ku.dk> Date: Thu, 15 Apr 2004 18:06:08 +0200 (CEST) X-Spam-Status: No, hits=-1.5 required=5.0 tests=BAYES_20,NO_REAL_NAME version=2.55 X-Spam-Level: X-Spam-Checker-Version: SpamAssassin 2.55 (1.174.2.19-2003-05-19-exp) Full_Name: Morten Welinder Version: snapshot OS: Submission from: (NULL) (65.213.85.218) Time to kick phyper's tires... The current version has very serious cancellation issues. For example, if you ask for a small right-tail you are likely to get total cancellation. For example phyper(59, 150, 150, 60, FALSE, FALSE) gives 6.372680161e-14. The right answer is dhyper(0, 150, 150, 60, FALSE) which is 5.111204798e-22. phyper is also really slow for large arguments. Therefore, I suggest using the code below. This is a sniplet from Gnumeric, so you'll have to s/gnm_float/double/ and s/gnum//, etc. The code isn't perfect. In fact, if i*(NR+NB) is close to n*NR, then this code can take a while. Not longer than the old code, though. /* * New phyper implementation. Copyright 2004 Morten Welinder. * Distributed under the GNU General Public License. * * Thanks to Ian Smith for ideas. */ /* * Calculate * * phyper (i, NR, NB, n, TRUE, FALSE) * [log] ---------------------------------- * dhyper (i, NR, NB, n, FALSE) * * without actually calling phyper. This assumes that * * i * (NR + NB) <= n * NR * */ static gnm_float pdhyper (gnm_float i, gnm_float NR, gnm_float NB, gnm_float n, gboolean log_p) { gnm_float sum = 0; gnm_float term = 1; while (i > 0 && term >= GNUM_EPSILON * sum) { term *= i * (NB - n + i) / (n + 1 - i) / (NR + 1 - i); sum += term; i--; } return log_p ? log1pgnum (sum) : 1 + sum; } gnm_float phyper (gnm_float i, gnm_float NR, gnm_float NB, gnm_float n, int lower_tail, int log_p) { gnm_float d, pd; #ifdef IEEE_754 if (isnangnum (i) || isnangnum (NR) || isnangnum (NB) || isnangnum (n)) return i + NR + NB + n; #endif i = floorgnum (i + 1e-7); NR = floorgnum (NR + 0.5); NB = floorgnum (NB + 0.5); n = floorgnum (n + 0.5); if (NR < 0 || NB < 0 || !finitegnum (NR + NB) || n < 0 || n > NR + NB) ML_ERR_return_NAN; if (i * (NR + NB) > n * NR) { /* Swap tails. */ gnm_float oldNB = NB; NB = NR; NR = oldNB; i = n - i - 1; lower_tail = !lower_tail; } if (i < 0) return R_DT_0; d = dhyper (i, NR, NB, n, log_p); pd = pdhyper (i, NR, NB, n, log_p); return log_p ? R_DT_log (d + pd) : R_D_Lval (d * pd); } ==> 6772.audit <== Fri Apr 16 18:26:36 2004 ripley moved from incoming to Accuracy ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Directory: Add-ons * PR# 974 * Subject: Lattice: panel.superpose with ordered factor groups From: John Maindonald Date: Sat, 9 Jun 2001 11:08:51 +1000 (EST) --The warning is standard S and R behaviour. --Probably xyplot needs to avoid it (by unclassing?) --Still there in lattice 0.3-0, 0.6-8 ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 974 <== From postmaster@franz.stat.wisc.edu Sat Jun 9 03:09:10 2001 Received: from franz.stat.wisc.edu (root@franz.stat.wisc.edu [128.105.174.95]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id DAA20814 for ; Sat, 9 Jun 2001 03:09:09 +0200 (MET DST) Received: from leonard.anu.edu.au (really [150.203.2.15]) by franz.stat.wisc.edu via smail with esmtp id (Debian Smail3.2.0.111) for ; Fri, 8 Jun 2001 20:09:01 -0500 (CDT) Received: (from u9801539@localhost) by leonard.anu.edu.au (8.9.3/8.9.3) id LAA01081; Sat, 9 Jun 2001 11:08:51 +1000 (EST) Date: Sat, 9 Jun 2001 11:08:51 +1000 (EST) From: John Maindonald Message-Id: <200106090108.LAA01081@leonard.anu.edu.au> To: r-bugs@r-project.org Subject: Lattice: panel.superpose with ordered factor groups Cc: deepayan@stat.wisc.edu # r-bugs@r-project.org The panel.superpose() function for xyplot() complains when the groups parameter is set to be an ordered factor. > library(mass) > data(cabbages) > library(lattice) > xyplot(HeadWt~VitC|Cult, panel=panel.superpose, + groups=Date, data=cabbages) > cabbages$Date <- ordered(cabbages$Date) # Make Date an ordered factor > xyplot(HeadWt~VitC|Cult, panel=panel.superpose, + groups=Date, data=cabbages) Warning messages: 1: Incompatible methods ("Ops.ordered", "Ops.factor") for "==" 2: Incompatible methods ("Ops.ordered", "Ops.factor") for "==" 3: Incompatible methods ("Ops.ordered", "Ops.factor") for "==" 4: Incompatible methods ("Ops.ordered", "Ops.factor") for "==" 5: Incompatible methods ("Ops.ordered", "Ops.factor") for "==" 6: Incompatible methods ("Ops.ordered", "Ops.factor") for "==" --please do not edit the information below-- Version: platform = i386-pc-mingw32 arch = x86 os = Win32 system = x86, Win32 status = major = 1 minor = 2.3 year = 2001 month = 04 day = 26 language = R Windows ME 4.90 (build 3000) Search Path: .GlobalEnv, package:grid, package:lattice, package:mass, package:ctest, Autoloads, package:base John Maindonald email: john.maindonald@anu.edu.au phone : (6125)3998 fax : (6125)5549 Statistical Consulting Unit, Room 1194, John Dedman Mathematical Sciences Building (Centre for Mathematics & Its Applications.) ==> 974.audit <== Sun Jun 17 09:39:36 2001 ripley changed notes Sun Jun 17 09:39:36 2001 ripley foobar Sun Jun 17 09:39:36 2001 ripley moved from incoming to Add-ons Sun Sep 30 15:13:56 2001 ripley changed notes Sun Sep 30 15:13:56 2001 ripley foobar Wed Jan 8 11:14:52 2003 ripley changed notes Wed Jan 8 11:14:52 2003 ripley foobar ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 1361 * Subject: Matrix identification bug From: hyu@stats.uwo.ca Date: Tue, 5 Mar 2002 21:19:46 +0100 (MET) --seems to be about Matrix package, not solve ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 1361 <== From hyu@stats.uwo.ca Tue Mar 5 21:19:46 2002 Received: from blueberry.kubism.ku.dk (blueberry [192.38.19.4]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id VAA08030 for ; Tue, 5 Mar 2002 21:19:46 +0100 (MET) Date: Tue, 5 Mar 2002 21:19:46 +0100 (MET) From: hyu@stats.uwo.ca Message-Id: <200203052019.VAA08030@pubhealth.ku.dk> To: R-bugs@biostat.ku.dk Subject: Matrix identification bug Full_Name: Hao Yu Version: 1.4.1 OS: Windws and Linux Submission from: (NULL) (129.100.45.161) Hi, Recently we did some benchmarks regarding matrix inverse operation and found some strange things. See the following results (the package Matrix is most updated). 1. load the function Toeplitz and library(Matrix) "Toeplitz" function(x) { matrix(x[1 + abs(outer(seq(along = x), seq(along = x), FUN = "-"))], byrow = T, ncol = length(x)) } 2. On a PIII 1.125GHz (256kb cache)desktop PC with win2k, the results are > system.time(solve(Toeplitz(.5^seq(0,999)))) [1] 38.25 0.22 42.01 NA NA > system.time(solve.Matrix(Toeplitz(.5^seq(0,999)))) [1] 28.01 0.17 30.94 NA NA On a PIII-M 1.2GHz (512kb cache) laptop PC with winxp pro, the results are > system.time(solve(Toeplitz(.5^seq(0,999)))) [1] 33.94 0.25 34.01 NA NA > system.time(solve.Matrix(Toeplitz(.5^seq(0,999)))) [1] 9.40 0.13 9.54 NA NA On a PIII Xeon 1GHz (256kb cache) RedHat 7.2 Linux, the results are > system.time(solve(Toeplitz(.5^seq(0,999)))) [1] 53.79 0.48 00.01 NA NA > system.time(solve.Matrix(Toeplitz(.5^seq(0,999)))) [1] 30.29 0.21 30.59 NA NA This was the first time that windows R outperformed linux R even PIII Xeon 1GHz is slightly faster than PIII 1.125GHz. Notice the abnormal result 9.4s from PIII-M 1.2GHz (cannot be three times faster). After using native fortran lapack (call to dpotrf and dpotri returns 9s for PIII Xeon), I realized that R may treat the matrix as a general matrix rather than symmetric except the abnormal case. To confirm this, I modified the R_LapackPP.cc in Matrix package and it turned out that is.MMatrix() failed to return true. Hence if (checkClass(classes, "Hermitian")) return new LaSymmMatDouble(x); is not executed in static LaMatDouble* asLaMatrix(SEXP x). Instead if (isMatrix(x)) { return new LaGenMatDouble(x); } is used. After replaced return new LaGenMatDouble(x); with return new LaSymmMatDouble(x); the result was (PIII Xeon 1GHz RedHat 7.2 Linux)) > system.time(solve.Matrix(Toeplitz(.5^seq(0,999)))) [1] 9.76 0.22 9.97 0.00 0.00 Clearly there are some bugs in solve() and solve.Matrix(). Notice that all have > is.Hermitian(Toeplitz(.5^seq(0,999))) [1] TRUE Thanks Hao Yu ==> 1361.followup.1 <== From ripley@stats.ox.ac.uk Tue Mar 5 22:32:30 2002 Received: from toucan.stats.ox.ac.uk (toucan.stats.ox.ac.uk [163.1.20.20]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id WAA08212 for ; Tue, 5 Mar 2002 22:32:30 +0100 (MET) Received: from tern.stats (pc2-oxfo4-0-cust120.oxf.cable.ntl.com [62.254.131.120]) by toucan.stats.ox.ac.uk (8.9.0/8.9.0) with ESMTP id VAA24644; Tue, 5 Mar 2002 21:31:57 GMT Date: Tue, 5 Mar 2002 21:32:31 +0000 (GMT Standard Time) From: Prof Brian D Ripley To: hyu@stats.uwo.ca cc: R-bugs@biostat.ku.dk Subject: Re: [Rd] Matrix identification bug (PR#1361) In-Reply-To: <200203052019.VAA08035@pubhealth.ku.dk> Message-ID: Sender: ripley@stats.ox.ac.uk MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII On Tue, 5 Mar 2002 hyu@stats.uwo.ca wrote: > Full_Name: Hao Yu > Version: 1.4.1 > OS: Windws and Linux > Submission from: (NULL) (129.100.45.161) [...] > the result was (PIII Xeon 1GHz RedHat 7.2 Linux)) > > system.time(solve.Matrix(Toeplitz(.5^seq(0,999)))) > [1] 9.76 0.22 9.97 0.00 0.00 > > Clearly there are some bugs in solve() and solve.Matrix(). Notice that all have Nothing in this report suggests a bug in solve(). Can you please provide some evidence for your claim (which is not `clear' to me), and > solve function (a, b, ...) UseMethod("solve") seems hard to get wrong! (It matters: solve is part of R: solve.Matrix is a contributed addon.) Incidentally, these sort of timings are crucially dependent on the BLAS library used, and in my experience that is a far larger factor that the OS used. However, that was not part of the report .... -- Brian D. Ripley, ripley@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 272860 (secr) Oxford OX1 3TG, UK Fax: +44 1865 272595 ==> 1361.audit <== Thu Mar 7 12:41:18 2002 ripley changed notes Thu Mar 7 12:41:18 2002 ripley foobar Thu Mar 7 12:41:18 2002 ripley moved from incoming to Add-ons ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 1662 * Subject: fisher.test FEXACT memory bug "should not occur" From: Martin Maechler Date: Thu, 13 Jun 2002 08:21:50 +0200 --The supplementary (table of sum one) is fixed for 1.5.1. --Detection code for the first problem has been added to 1.5.1 which will stop --the crash, but the underlying cause is still open. ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 1662 <== From maechler@stat.math.ethz.ch Thu Jun 13 08:22:43 2002 Received: from stat.math.ethz.ch (daemon@hypatia.ethz.ch [129.132.58.23]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id IAA16033 for ; Thu, 13 Jun 2002 08:22:42 +0200 (MET DST) Received: (from daemon@localhost) by stat.math.ethz.ch (8.9.1/8.9.1) id IAA01485 for ; Thu, 13 Jun 2002 08:21:35 +0200 (MET DST) Received: from lynne(129.132.58.30), claiming to be "lynne.ethz.ch" via SMTP by hypatia, id smtpdAAAa000Jg; Thu Jun 13 08:21:26 2002 Received: (from maechler@localhost) by lynne.ethz.ch (8.11.6/8.11.2) id g5D6LoC13482; Thu, 13 Jun 2002 08:21:50 +0200 X-Authentication-Warning: lynne.ethz.ch: maechler set sender to maechler@stat.math.ethz.ch using -f From: Martin Maechler MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Transfer-Encoding: 7bit Message-ID: <15624.14846.214663.363944@gargle.gargle.HOWL> Date: Thu, 13 Jun 2002 08:21:50 +0200 To: R-bugs@stat.math.ethz.ch Subject: fisher.test FEXACT memory bug "should not occur" X-Mailer: VM 7.07 under Emacs 20.7.1 Reply-To: Martin Maechler This is a bad bug as reported by Robin Hankin, it is still in "R-patched" ... ##- From: Robin Hankin ##- To: r-help@stat.math.ethz.ch ##- Subject: [R] possum sleeping: thanks and fisher.test() FEXACT error ##- Date: Thu, 13 Jun 2002 16:46:26 +1200 ## ..... ## Example slighlty modified (MM) d4 <- matrix(c(0, 0, 0, 0, 0, 0, 3, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 2, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 4, 0, 2, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 2, 0, 0, 0, 2, 2, 0, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 2, 0, 0, 0, 3, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 2, 0, 0, 0, 0, 8, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 1, 0, 0, 2, 0, 0, 0, 0, 2, 0, 0, 1, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 4, 0, 0, 0, 0), nr=50) fisher.test(d4) ##- Error in fisher.test(alldata[, 2:5]) : FEXACT error 30. ##- Stack length exceeded in f3xact. ##- This problem should not occur. ==> 1662.followup.1 <== From ligges@amadeus.statistik.uni-dortmund.de Thu Jun 13 08:58:58 2002 Received: from nx5.HRZ.Uni-Dortmund.DE (nx5.HRZ.Uni-Dortmund.DE [129.217.131.21]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id IAA16332 for ; Thu, 13 Jun 2002 08:58:57 +0200 (MET DST) Received: from amadeus.statistik.uni-dortmund.de by nx5.HRZ.Uni-Dortmund.DE via smtp-local with ESMTP; Thu, 13 Jun 2002 08:56:18 +0200 Received: from statistik.uni-dortmund.de ([129.217.207.201]) by amadeus.statistik.uni-dortmund.de (8.11.6+Sun/8.11.6) with ESMTP id g5D6uH708478; Thu, 13 Jun 2002 08:56:18 +0200 (MET DST) Message-ID: <3D084237.3F455ED2@statistik.uni-dortmund.de> Date: Thu, 13 Jun 2002 08:56:55 +0200 From: Uwe Ligges Organization: Fachbereich Statistik, Universitaet Dortmund X-Mailer: Mozilla 4.79 [en] (WinNT; U) X-Accept-Language: en MIME-Version: 1.0 To: maechler@stat.math.ethz.ch CC: R-bugs@biostat.ku.dk Subject: Re: fisher.test FEXACT memory bug "should not occur" (PR#1662) References: <200206130622.IAA16038@pubhealth.ku.dk> Content-Type: text/plain; charset=us-ascii Content-Transfer-Encoding: 7bit maechler@stat.math.ethz.ch wrote: > > This is a bad bug as reported by Robin Hankin, > it is still in "R-patched" ... > > ##- From: Robin Hankin > ##- To: r-help@stat.math.ethz.ch > ##- Subject: [R] possum sleeping: thanks and fisher.test() FEXACT error > ##- Date: Thu, 13 Jun 2002 16:46:26 +1200 > > ## ..... > > ## Example slighlty modified (MM) > > d4 <- matrix(c(0, 0, 0, 0, 0, 0, 3, 0, 1, 0, 0, 0, 0, 0, 0, > 1, 0, 0, 0, 0, 1, 0, 0, 2, 0, 0, 0, 1, 0, 0, > 0, 1, 0, 1, 0, 4, 0, 2, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, > 0, 1, 0, 0, 2, 0, 0, 0, 2, 2, 0, 1, 0, 0, 0, > 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, > 1, 0, 0, 0, 2, 0, 0, 0, 3, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, > 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, > 0, 0, 0, 0, 2, 0, 0, 0, 0, 8, 0, 0, 0, 3, 0, > 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, > 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 1, 0, > 0, 2, 0, 0, 0, 0, 2, 0, 0, 1, 3, 0, 0, 0, 0, > 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 4, 0, 0, 0, 0), > nr=50) > > fisher.test(d4) > ##- Error in fisher.test(alldata[, 2:5]) : FEXACT error 30. > ##- Stack length exceeded in f3xact. > ##- This problem should not occur. Just for the record: repeat try(fisher.test(d4)) results in a crash after a few iterations (R-1.5.0 patched on WinNT 4.0). Uwe Ligges ==> 1662.reply.1 <== From maechler@stat.math.ethz.ch Thu Jun 13 12:17:02 2002 Received: from stat.math.ethz.ch (daemon@hypatia.ethz.ch [129.132.58.23]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id MAA19850 for ; Thu, 13 Jun 2002 12:17:01 +0200 (MET DST) Received: (from daemon@localhost) by stat.math.ethz.ch (8.9.1/8.9.1) id MAA05039; Thu, 13 Jun 2002 12:15:37 +0200 (MET DST) Received: from lynne(129.132.58.30), claiming to be "lynne.ethz.ch" via SMTP by hypatia, id smtpdAAAa001Eg; Thu Jun 13 12:15:28 2002 Received: (from maechler@localhost) by lynne.ethz.ch (8.11.6/8.11.2) id g5DAFsK19920; Thu, 13 Jun 2002 12:15:54 +0200 X-Authentication-Warning: lynne.ethz.ch: maechler set sender to maechler@stat.math.ethz.ch using -f From: Martin Maechler MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Transfer-Encoding: 7bit Message-ID: <15624.28889.858685.234909@gargle.gargle.HOWL> Date: Thu, 13 Jun 2002 12:15:53 +0200 To: Prof Brian D Ripley Cc: R-bugs@stat.math.ethz.ch, R-core@stat.math.ethz.ch Subject: Re: bad fisher.test() bug (PR#1662) In-Reply-To: References: <15624.23861.930057.723845@gargle.gargle.HOWL> X-Mailer: VM 7.07 under Emacs 20.7.1 Reply-To: Martin Maechler (CC'ed to R-bugs ``for the record'') >>>>> "BDR" == Prof Brian D Ripley writes: BDR> On Thu, 13 Jun 2002, Martin Maechler wrote: >> >>>>> "MM" == Martin Maechler >> writes: >> >> >>>>> "BDR" == Brian D Ripley >> writes: BDR> Martin, What makes this a `bad bug'? Are you getting a BDR> seg fault? >> MM> yes, always in the half a dozen restarts I tried. >> (in the mean time I had one case where it did not ..) >> BDR> Like Uwe this works for me the first few times and I do BDR> then get a backtrace (in memory.c, so this is almost BDR> certainly an earlier overrun). >> MM> I see. This has been different for me. >> BDR> I think we've had problems before with FEXACT BDR> incorrectly specifying the required sizes of its BDR> workspaces. >> (I'm not so sure this is the problem here.) >> >> I found that >> >> fisher.test(d4[1:30,]) >> >> gives a direct segmentation fault (without any error >> message) for me. >> BDR> I've had to compile up gdb 5.2 to get the correct BDR> information out of R on Linux compiled under gcc-3.1. BDR> That shows (on the original problem) BDR> Loaded symbols for BDR> /users/ripley/R/R-patched/library/ctest/libs/ctest.so BDR> #0 f3xact (nrow=0x40051074, irow=0x4005113c, BDR> ncol=0x40051000, icol=0x40051120, dlp=0x4041b110, BDR> mm=0x40050fc0, fact=0x4035f020, ico=0x4035f4dc, BDR> iro=0x4035f5a4, it=0x4035f66c, lb=0x4035f734, BDR> nr=0x4035f7fc, nt=0x4035f8c4, nu=0x4035f98c, BDR> itc=0x4035fa54, ist=0x40360094, stv=0x40364f18, BDR> alen=0x40365ba0, tol=0x4004f878) at BDR> /users/ripley/R/cvs/R-patched/src/library/ctest/src/fexact.c:1125 BDR> 1125 ist[itp] = -1; (gdb) print itp $1 = -3544 BDR> The same with your variant (except the value of itp). Okay, I found a much smaller table --- and a different place for the seg.fault : > fisher.test(cbind(0,c(0,0,1))) Program received signal SIGSEGV, Segmentation fault. f2xact (nrow=0x8a12ec0, ncol=0x8a12ee0, table=0x8b7bb38, ldtabl=0x8a12f00, expect=0x89dad68, percnt=0x89dad90, emin=0x89dadb8, prt=0x89dade0, pre=0x89dae08, fact=0x403b3020, ico=0x403b302c, iro=0x403b3038, kyy=0x403b3044, idif=0x403b3050, irn=0x403b3058, key=0x403b49d4, ldkey=0xbfffc398, ipoin=0x403b5d44, stp=0x403b70b0, ldstp=0xbfffc39c, ifrq=0x403ffef4, dlp=0x4046d450, dsp=0x4046fb30, tm=0x40472210, key2=0x404748f4, iwk=0x403b3060, rwk=0x403b3d30) at ../../../../../R/src/library/ctest/src/fexact.c:524 524 obs += fact[ico[j]] - dd; Martin ==> 1662.followup.2 <== From ripley@stats.ox.ac.uk Thu Jun 13 13:16:27 2002 Received: from toucan.stats.ox.ac.uk (toucan.stats.ox.ac.uk [163.1.20.20]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id NAA20855 for ; Thu, 13 Jun 2002 13:16:27 +0200 (MET DST) Received: from auk.stats (pc1-oxfd1-4-cust16.oxf.cable.ntl.com [62.254.131.16]) by toucan.stats.ox.ac.uk (8.9.0/8.9.0) with ESMTP id MAA15227; Thu, 13 Jun 2002 12:15:44 +0100 (BST) Date: Thu, 13 Jun 2002 12:14:46 +0100 (BST) From: Prof Brian D Ripley Sender: ripley@auk.stats To: maechler@stat.math.ethz.ch cc: R-bugs@biostat.ku.dk, Subject: Re: fisher.test FEXACT memory bug "should not occur" (PR#1662) In-Reply-To: <200206130622.IAA16038@pubhealth.ku.dk> Message-ID: MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII This occurred because integer overflow was giving negative keys. Append the following line at line 1083 of fexact.c if (ipn < 1) ipn += ldst; /* because key might be negative */ > fisher.test(d4) Fisher's Exact Test for Count Data data: d4 p-value = < 2.2e-16 alternative hypothesis: two.sided On Thu, 13 Jun 2002 maechler@stat.math.ethz.ch wrote: > This is a bad bug as reported by Robin Hankin, > it is still in "R-patched" ... > > ##- From: Robin Hankin > ##- To: r-help@stat.math.ethz.ch > ##- Subject: [R] possum sleeping: thanks and fisher.test() FEXACT error > ##- Date: Thu, 13 Jun 2002 16:46:26 +1200 > > ## ..... > > ## Example slighlty modified (MM) > > d4 <- matrix(c(0, 0, 0, 0, 0, 0, 3, 0, 1, 0, 0, 0, 0, 0, 0, > 1, 0, 0, 0, 0, 1, 0, 0, 2, 0, 0, 0, 1, 0, 0, > 0, 1, 0, 1, 0, 4, 0, 2, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, > 0, 1, 0, 0, 2, 0, 0, 0, 2, 2, 0, 1, 0, 0, 0, > 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, > 1, 0, 0, 0, 2, 0, 0, 0, 3, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, > 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, > 0, 0, 0, 0, 2, 0, 0, 0, 0, 8, 0, 0, 0, 3, 0, > 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, > 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 1, 0, > 0, 2, 0, 0, 0, 0, 2, 0, 0, 1, 3, 0, 0, 0, 0, > 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 4, 0, 0, 0, 0), > nr=50) > > fisher.test(d4) > ##- Error in fisher.test(alldata[, 2:5]) : FEXACT error 30. > ##- Stack length exceeded in f3xact. > ##- This problem should not occur. > > > -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- > r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html > Send "info", "help", or "[un]subscribe" > (in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch > _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ > -- Brian D. Ripley, ripley@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 272860 (secr) Oxford OX1 3TG, UK Fax: +44 1865 272595 ==> 1662.followup.3 <== From ripley@stats.ox.ac.uk Thu Jun 13 17:57:14 2002 Received: from stat.math.ethz.ch (daemon@hypatia.ethz.ch [129.132.58.23]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id RAA24468 for ; Thu, 13 Jun 2002 17:57:14 +0200 (MET DST) From: ripley@stats.ox.ac.uk Received: (from daemon@localhost) by stat.math.ethz.ch (8.9.1/8.9.1) id RAA20365; Thu, 13 Jun 2002 17:56:07 +0200 (MET DST) Received: from toucan.stats.ox.ac.uk(163.1.20.20) via SMTP by hypatia, id smtpdAAAa004wF; Thu Jun 13 17:55:59 2002 Received: from gannet.stats (gannet.stats [163.1.20.127]) by toucan.stats.ox.ac.uk (8.9.0/8.9.0) with ESMTP id QAA16075; Thu, 13 Jun 2002 16:56:21 +0100 (BST) Date: Thu, 13 Jun 2002 16:56:20 +0100 (BST) X-X-Sender: To: Martin Maechler cc: Subject: Re: bad fisher.test() bug (PR#1662) In-Reply-To: <15624.28889.858685.234909@gargle.gargle.HOWL> Message-ID: MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII On Thu, 13 Jun 2002, Martin Maechler wrote: > Okay, I found a much smaller table --- and a different place > for the seg.fault : > > fisher.test(cbind(0,c(0,0,1))) That one is easy: You (MM) have calculated fact[2], but fact is only of length 2 (for values 0 and 1), so ico gets overwritten. I now suspect that my fix for the original bug is merely a fix for a symptom, so am going to replace it by an error message. -- Brian D. Ripley, ripley@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 272860 (secr) Oxford OX1 3TG, UK Fax: +44 1865 272595 ==> 1662.audit <== Thu Jun 13 19:15:42 2002 ripley changed notes Thu Jun 13 19:15:42 2002 ripley foobar Thu Jun 13 19:15:42 2002 ripley moved from incoming to Add-ons ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 1729 * Subject: problem with qq( ) From: Jarno.Tuimala@Helsinki.Fi Date: Tue, 2 Jul 2002 10:37:10 +0200 (MET DST) --A report on lattice ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 1729 <== From Jarno.Tuimala@Helsinki.Fi Tue Jul 2 10:37:11 2002 Received: from blueberry (blueberry [192.38.19.4]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id KAA28248 for ; Tue, 2 Jul 2002 10:37:10 +0200 (MET DST) Date: Tue, 2 Jul 2002 10:37:10 +0200 (MET DST) From: Jarno.Tuimala@Helsinki.Fi Message-Id: <200207020837.KAA28248@pubhealth.ku.dk> To: R-bugs@biostat.ku.dk Subject: problem with qq( ) Full_Name: Jarno Tuimala Version: 1.5.0 OS: Windows Nt Submission from: (NULL) (193.166.1.21) Running the following analysis gives as a result kukot vs. kisut. aikaero <- -27.5+5*(0:10) frekvenssi.m <- c(0,1,15,37,53,45,23,18,4,1,1) frekvenssi.n <- c(2,8,12,19,33,47,42,15,2,1,0) win.graph() aikaero <- rep(c(aikaero,aikaero), c(frekvenssi.m,frekvenssi.n)) a.ero <- data.frame(aikaero, sukupuoli=rep(c("kukot","kisut"), c(sum(frekvenssi.m),sum(frekvenssi.n)))) qq(sukupuoli ~ aikaero, data=a.ero) However, running the following analysis, gives miehet vs. miehet. aikaero <- -27.5+5*(0:10) frekvenssi.m <- c(0,1,15,37,53,45,23,18,4,1,1) frekvenssi.n <- c(2,8,12,19,33,47,42,15,2,1,0) win.graph() aikaero <- rep(c(aikaero,aikaero), c(frekvenssi.m,frekvenssi.n)) a.ero <- data.frame(aikaero, sukupuoli=rep(c("miehet","naiset"), c(sum(frekvenssi.m),sum(frekvenssi.n)))) qq(sukupuoli ~ aikaero, data=a.ero) This seems to happen everytime "sukupuoli" contains one variable name, which starts with a letter "n". If both names start with "n" or other letters the results is miehet vs naiset as expected. So this works fine: a.ero <- data.frame(aikaero, sukupuoli=rep(c("nisset","naiset"), c(sum(frekvenssi.m),sum(frekvenssi.n)))) qq(sukupuoli ~ aikaero, data=a.ero) And this also: a.ero <- data.frame(aikaero, sukupuoli=rep(c("miehet","maiset"), c(sum(frekvenssi.m),sum(frekvenssi.n)))) qq(sukupuoli ~ aikaero, data=a.ero) Where's the problem? ==> 1729.followup.1 <== From ripley@stats.ox.ac.uk Tue Jul 2 11:19:02 2002 Received: from toucan.stats.ox.ac.uk (toucan.stats.ox.ac.uk [163.1.20.20]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id LAA28796 for ; Tue, 2 Jul 2002 11:19:02 +0200 (MET DST) From: ripley@stats.ox.ac.uk Received: from gannet.stats (gannet.stats [163.1.20.127]) by toucan.stats.ox.ac.uk (8.9.0/8.9.0) with ESMTP id KAA21799; Tue, 2 Jul 2002 10:18:19 +0100 (BST) Date: Tue, 2 Jul 2002 10:18:18 +0100 (BST) X-X-Sender: To: cc: , Subject: Re: problem with qq( ) (PR#1729) In-Reply-To: <200207020837.KAA28253@pubhealth.ku.dk> Message-ID: MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII 1) This is a report on the lattice package, and we both need you to mention that and to give the version number of lattice that you used to be sure of beeing able to reproduce this. 2) It seems that the plots are actually correct, just the labels are wrong. You can correct that by specifying xlab and ylab in the call. 3) The error appears to be a typo: if (missing(ylab)) ylab <- if (is.f.y) unique(levels(y))[y] else paste("y:", as.character(unique(levels(y)[[2]]))) should be [2] for consistency. On Tue, 2 Jul 2002 Jarno.Tuimala@Helsinki.Fi wrote: > Full_Name: Jarno Tuimala > Version: 1.5.0 > OS: Windows Nt > Submission from: (NULL) (193.166.1.21) > > > Running the following analysis gives as a result kukot vs. kisut. > > aikaero <- -27.5+5*(0:10) > frekvenssi.m <- c(0,1,15,37,53,45,23,18,4,1,1) > frekvenssi.n <- c(2,8,12,19,33,47,42,15,2,1,0) > win.graph() > aikaero <- rep(c(aikaero,aikaero), c(frekvenssi.m,frekvenssi.n)) > a.ero <- data.frame(aikaero, > sukupuoli=rep(c("kukot","kisut"), > c(sum(frekvenssi.m),sum(frekvenssi.n)))) > qq(sukupuoli ~ aikaero, data=a.ero) > > However, running the following analysis, gives miehet vs. miehet. > > aikaero <- -27.5+5*(0:10) > frekvenssi.m <- c(0,1,15,37,53,45,23,18,4,1,1) > frekvenssi.n <- c(2,8,12,19,33,47,42,15,2,1,0) > win.graph() > aikaero <- rep(c(aikaero,aikaero), c(frekvenssi.m,frekvenssi.n)) > a.ero <- data.frame(aikaero, > sukupuoli=rep(c("miehet","naiset"), > c(sum(frekvenssi.m),sum(frekvenssi.n)))) > qq(sukupuoli ~ aikaero, data=a.ero) > > This seems to happen everytime "sukupuoli" contains one variable name, which > starts with a letter "n". If both names start with "n" or other letters the > results is miehet vs naiset as expected. > > So this works fine: > > a.ero <- data.frame(aikaero, > sukupuoli=rep(c("nisset","naiset"), > c(sum(frekvenssi.m),sum(frekvenssi.n)))) > qq(sukupuoli ~ aikaero, data=a.ero) > > And this also: > > a.ero <- data.frame(aikaero, > sukupuoli=rep(c("miehet","maiset"), > c(sum(frekvenssi.m),sum(frekvenssi.n)))) > qq(sukupuoli ~ aikaero, data=a.ero) > > Where's the problem? > > > -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- > r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html > Send "info", "help", or "[un]subscribe" > (in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch > _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ > -- Brian D. Ripley, ripley@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 272860 (secr) Oxford OX1 3TG, UK Fax: +44 1865 272595 ==> 1729.followup.2 <== From deepayansarkar@yahoo.com Mon Jul 8 19:44:52 2002 Received: from web13902.mail.yahoo.com (web13902.mail.yahoo.com [216.136.175.28]) by pubhealth.ku.dk (8.9.3/8.9.1) with SMTP id TAA03638 for ; Mon, 8 Jul 2002 19:44:51 +0200 (MET DST) Message-ID: <20020708174408.63945.qmail@web13902.mail.yahoo.com> Received: from [202.54.53.232] by web13902.mail.yahoo.com via HTTP; Mon, 08 Jul 2002 10:44:08 PDT Date: Mon, 8 Jul 2002 10:44:08 -0700 (PDT) From: Deepayan Sarkar Subject: Re: problem with qq( ) (PR#1729) To: R-bugs@biostat.ku.dk In-Reply-To: MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii --- ripley@stats.ox.ac.uk wrote: > 1) This is a report on the lattice package, and we both need you to > mention that and to give the version number of lattice that you used > to be sure of beeing able to reproduce this. > > 2) It seems that the plots are actually correct, just the labels are > wrong. You can correct that by specifying xlab and ylab in the call. > > 3) The error appears to be a typo: > > if (missing(ylab)) > ylab <- if (is.f.y) > unique(levels(y))[y] > else paste("y:", as.character(unique(levels(y)[[2]]))) > > should be [2] for consistency. Thanks, fixed for future releases. > On Tue, 2 Jul 2002 Jarno.Tuimala@Helsinki.Fi wrote: > > > Full_Name: Jarno Tuimala > > Version: 1.5.0 > > OS: Windows Nt > > Submission from: (NULL) (193.166.1.21) > > > > > > Running the following analysis gives as a result kukot vs. kisut. > > > > aikaero <- -27.5+5*(0:10) > > frekvenssi.m <- c(0,1,15,37,53,45,23,18,4,1,1) > > frekvenssi.n <- c(2,8,12,19,33,47,42,15,2,1,0) > > win.graph() > > aikaero <- rep(c(aikaero,aikaero), c(frekvenssi.m,frekvenssi.n)) > > a.ero <- data.frame(aikaero, > > sukupuoli=rep(c("kukot","kisut"), > > c(sum(frekvenssi.m),sum(frekvenssi.n)))) > > qq(sukupuoli ~ aikaero, data=a.ero) > > > > However, running the following analysis, gives miehet vs. miehet. > > > > aikaero <- -27.5+5*(0:10) > > frekvenssi.m <- c(0,1,15,37,53,45,23,18,4,1,1) > > frekvenssi.n <- c(2,8,12,19,33,47,42,15,2,1,0) > > win.graph() > > aikaero <- rep(c(aikaero,aikaero), c(frekvenssi.m,frekvenssi.n)) > > a.ero <- data.frame(aikaero, > > sukupuoli=rep(c("miehet","naiset"), > > c(sum(frekvenssi.m),sum(frekvenssi.n)))) > > qq(sukupuoli ~ aikaero, data=a.ero) > > > > This seems to happen everytime "sukupuoli" contains one variable name, > which > > starts with a letter "n". If both names start with "n" or other letters the > > results is miehet vs naiset as expected. > > > > So this works fine: > > > > a.ero <- data.frame(aikaero, > > sukupuoli=rep(c("nisset","naiset"), > > c(sum(frekvenssi.m),sum(frekvenssi.n)))) > > qq(sukupuoli ~ aikaero, data=a.ero) > > > > And this also: > > > > a.ero <- data.frame(aikaero, > > sukupuoli=rep(c("miehet","maiset"), > > c(sum(frekvenssi.m),sum(frekvenssi.n)))) > > qq(sukupuoli ~ aikaero, data=a.ero) > > > > Where's the problem? > > > > > > > -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- > > r-devel mailing list -- Read > http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html > > Send "info", "help", or "[un]subscribe" > > (in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch > > > _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ > > > > -- > Brian D. Ripley, ripley@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 272860 (secr) > Oxford OX1 3TG, UK Fax: +44 1865 272595 > > -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- > r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html > Send "info", "help", or "[un]subscribe" > (in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch > _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._ __________________________________________________ Do You Yahoo!? Sign up for SBC Yahoo! Dial - First Month Free http://sbc.yahoo.com ==> 1729.audit <== Tue Jul 2 21:14:42 2002 ripley moved from incoming to Add-ons Sun Aug 25 11:58:22 2002 ripley changed notes Sun Aug 25 11:58:22 2002 ripley foobar ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 1972 * Subject: lattice install From: robert.king@newcastle.edu.au Date: Mon, 2 Sep 2002 04:55:41 +0200 (MET DST) --Perhaps lattice should require(grid) and print a clearer message? ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 1972 <== From robert.king@newcastle.edu.au Mon Sep 2 04:55:43 2002 Received: from blueberry (blueberry [192.38.19.4]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id EAA26744 for ; Mon, 2 Sep 2002 04:55:41 +0200 (MET DST) Date: Mon, 2 Sep 2002 04:55:41 +0200 (MET DST) From: robert.king@newcastle.edu.au Message-Id: <200209020255.EAA26744@pubhealth.ku.dk> To: R-bugs@biostat.ku.dk Subject: lattice install Full_Name: Robert King Version: 1.3 OS: windows Submission from: (NULL) (134.148.4.19) This looks like a problem with package dependencies. Should be sending this to the package author, not r-bugs? An attempt to install lattice from the windows binary runs OK, but attempting to load the library an error occurs: > install.packages("D:/Software/R/lattice.zip", .lib.loc[1], CRAN = NULL) updating HTML package descriptions > library(lattice) Error in parse(file, n, text, prompt) : syntax error on line 2022 I then tried installing under linux with R 1.5.1 This installed OK, but loading the library gave another error: > library(lattice) Error in library(grid) : There is no package called `grid' Installing grid on linux resulted in lattice loading and working properly Installing grid also fixed the problem on windows. ==> 1972.audit <== Tue Sep 3 19:47:36 2002 thomas changed notes Tue Sep 3 19:47:36 2002 thomas foobar Tue Sep 3 19:47:36 2002 thomas moved from incoming to Graphics Wed Jan 8 11:27:36 2003 ripley foobar Wed Jan 8 11:27:36 2003 ripley moved from Graphics to Add-ons ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 1974 * Subject: Rwave installation problem From: ld-temp-qt3i@pobox.com Date: Mon, 2 Sep 2002 09:27:36 +0200 (MET DST) --Instead, this should use ISO C headers, namely ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 1974 <== From ld-temp-qt3i@pobox.com Mon Sep 2 09:27:37 2002 Received: from blueberry (blueberry [192.38.19.4]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id JAA27886 for ; Mon, 2 Sep 2002 09:27:36 +0200 (MET DST) Date: Mon, 2 Sep 2002 09:27:36 +0200 (MET DST) From: ld-temp-qt3i@pobox.com Message-Id: <200209020727.JAA27886@pubhealth.ku.dk> To: R-bugs@biostat.ku.dk Subject: Rwave installation problem Full_Name: Linc Davis Version: 1.5.1 OS: Mac OS 10.1.5 Submission from: (NULL) (68.4.90.94) When running install.packages("Rwave") as root, I got this error: In file included from compinteg.c:12: kernel.h:22: malloc.h: No such file or directory make: *** [compinteg.o] Error 1 ERROR: compilation failed for package 'Rwave' There are 3 header files 'malloc.h' in Mac OS X: one for generic C, one for Objective-C, and a private one for the OS. The R installation script looks for C headers in /usr/include, but the actual path to this one is /usr/include/sys/malloc.h I worked around the problem by copying the file as indicated, then I deleted the copy after compiling the package. To fix the bug, the script should be changed to look for C headers in both /usr/include and /usr/include/sys. Thanks. ==> 1974.audit <== Tue Sep 17 11:20:29 2002 maechler foobar Tue Sep 17 11:20:29 2002 maechler moved from incoming to Add-ons Wed Nov 20 16:44:15 2002 ripley changed notes Wed Nov 20 16:44:15 2002 ripley foobar ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 2173 * Subject: xlim in plot.survfit() [with a discussion on "..."] From: jerome@hivnet.ubc.ca Date: Wed, 16 Oct 2002 18:46:11 +0200 (MET DST) -- no notes -- ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 2173 <== From jerome@hivnet.ubc.ca Wed Oct 16 18:46:12 2002 Received: from blueberry (blueberry [192.38.19.4]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id SAA00603 for ; Wed, 16 Oct 2002 18:46:11 +0200 (MET DST) Date: Wed, 16 Oct 2002 18:46:11 +0200 (MET DST) From: jerome@hivnet.ubc.ca Message-Id: <200210161646.SAA00603@pubhealth.ku.dk> To: R-bugs@biostat.ku.dk Subject: xlim in plot.survfit() [with a discussion on "..."] Full_Name: Jerome Asselin Version: 1.6.0 OS: RedHat 7.2 Submission from: (NULL) (24.83.203.63) Hello, I am trying to draw a legend on top of survival curves using plot.survfit(). As in the example below, I would like to specify a large interval for the x-axis. I can achieve such result using "xlim". However, an error occurs if I use the legend.pos and legend.text parameters as well. library(survival) TIME <- Surv(c(23,52,12,56,23,21),c(T,F,T,F,F,T)) par(mfrow=c(2,1)) plot(survfit(TIME~1),conf.int=F,xlim=c(0,100), main="xlim=c(0,100) without legend") plot(survfit(TIME~1),conf.int=F,xlim=c(0,100), legend.pos=0,legend.text="Survival Curve", main="xlim=c(0,100) with legend (ERROR!)") It looks like plot.survfit() attempts to send the "xlim" parameter to legend(), which does not accept "xlim" as a parameter. Hence, an error occurs and the legend is not plotted. I believe the problem is due to the fact that the additional parameters "..." that are supplied to plot.survfit() do not necessarily match the range of possible parameters for both plot() and legend(). In my case, "xlim" is matched by plot(), but not by legend(). I know how to plot the legend by directly using the legend() function. However, in order to correct the small problem with plot.survfit(), I would suggest that plot.survfit() be modified to perform an exact match of the parameters in "..." onto legend(). Perhaps there is already a function to perform such a task. If not, I believe it could provide a practical solution for dealing with "..." and subfunctions that do not have a built-in "..." parameter. That's just a suggestion, though. Thank you very much! Jerome Asselin ==> 2173.audit <== Wed Oct 16 22:58:03 2002 ripley moved from incoming to Add-ons ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 2322 * Subject: simplex From: george@lecompte.org Date: Sat, 23 Nov 2002 17:30:37 +0100 (MET) --report on boot, I think (not mentioned, though) ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 2322 <== From george@lecompte.org Sat Nov 23 17:30:38 2002 Received: from blueberry (blueberry [192.38.19.4]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id RAA16415 for ; Sat, 23 Nov 2002 17:30:37 +0100 (MET) Date: Sat, 23 Nov 2002 17:30:37 +0100 (MET) From: george@lecompte.org Message-Id: <200211231630.RAA16415@pubhealth.ku.dk> To: R-bugs@biostat.ku.dk Subject: simplex Full_Name: George LeCompte Version: 1.6.1 (windows) OS: windows98 Submission from: (NULL) (208.8.162.135) simplex won't work with only >= constraints examples follow: this works: simplex(a=1,A1=1,b1=10,A2=1,b2=1) but this doesn't simplex(1,A2=1,b2=1) Linear Programming Results Call : simplex(a = 1, A2 = 1, b2 = 1) Minimization Problem with Objective Function Coefficients x1 1 No feasible solution could be found ==> 2322.audit <== Tue Nov 26 21:45:55 2002 ripley changed notes Tue Nov 26 21:45:55 2002 ripley foobar Tue Nov 26 21:45:55 2002 ripley moved from incoming to Add-ons ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 2352 * Subject: avas: segmentation fault on empty args From: vograno@arbitrade.com Date: Fri, 6 Dec 2002 19:41:48 +0100 (MET) -- no notes -- ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 2352 <== From vograno@arbitrade.com Fri Dec 6 19:41:49 2002 Received: from blueberry (blueberry [192.38.19.4]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id TAA11396 for ; Fri, 6 Dec 2002 19:41:48 +0100 (MET) Date: Fri, 6 Dec 2002 19:41:48 +0100 (MET) From: vograno@arbitrade.com Message-Id: <200212061841.TAA11396@pubhealth.ku.dk> To: R-bugs@biostat.ku.dk Subject: avas: segmentation fault on empty args Full_Name: Vadim Ogranovich Version: 1.6.0 OS: Red Hat 7.1 Submission from: (NULL) (209.99.241.1) Segmentation fault occurs when empty vectors passed to avas as arguments. > library("acepack") library("acepack") > avas(numeric(0),numeric(0)) avas(numeric(0),numeric(0)) Process R segmentation fault at Fri Dec 6 10:31:06 2002 ==> 2352.audit <== Sat Dec 7 09:32:42 2002 ripley foobar Sat Dec 7 09:32:42 2002 ripley moved from incoming to Add-ons ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 2369 * Subject: convergence problem with nlme() From: jerome@hivnet.ubc.ca Date: Thu, 12 Dec 2002 23:43:14 +0100 (MET) -- no notes -- ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 2369 <== From jerome@hivnet.ubc.ca Thu Dec 12 23:43:15 2002 Received: from blueberry (blueberry [192.38.19.4]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id XAA27889 for ; Thu, 12 Dec 2002 23:43:14 +0100 (MET) Date: Thu, 12 Dec 2002 23:43:14 +0100 (MET) From: jerome@hivnet.ubc.ca Message-Id: <200212122243.XAA27889@pubhealth.ku.dk> To: R-bugs@biostat.ku.dk Subject: convergence problem with nlme() Full_Name: Jerome Asselin Version: 1.6.1 OS: linux redhat 7.2 Submission from: (NULL) (142.103.173.179) I am using the nlme package version 3.1-33. I tried an example from the file "library/nlme/scripts/ch08.R" that comes with the nlme package but I did not obtain convergence. I assume that this example worked at some point in the past, but I cannot determine why it is not working anymore. So I am unable to determine whether this is a numerical bug or a documentation bug (in the file "ch08.R"). More specifically, I tried the example below (in which I added the option msVerbose=T.) library(nlme) data(Quinidine) control <- nlmeControl(msVerbose = T) fm1Quin.nlme <- nlme(conc ~ quinModel(Subject, time, conc, dose, interval, lV, lKa, lCl), data = Quinidine, fixed = lV + lKa + lCl ~ 1, random = pdDiag(lV + lCl ~ 1), groups = ~ Subject, start = list(fixed = c(5, -0.3, 2)), na.action = NULL, naPattern = ~ !is.na(conc), control = control) #Error in nlme.formula(conc ~ quinModel(Subject, time, conc, dose, interval, : # Maximum number of iterations reached without convergence Below is the trace of the calculations by nlme(). iteration = 0 Step: [1] 0 0 Parameter: [1] 0.9645714 1.2393630 Function Value [1] 1053.026 Gradient: [1] -5.1732725 -0.3194242 iteration = 10 Parameter: [1] 5.492493 1.275746 Function Value [1] 1050.678 Gradient: [1] -0.0004125939 0.0012034419 Iteration limit exceeded. Algorithm failed. iteration = 0 Step: [1] 0 0 Parameter: [1] 5.492303 0.860876 Function Value [1] 1047.32 Gradient: [1] 8.080429e-04 -8.308982e-06 iteration = 6 Parameter: [1] -1772.4732661 0.3624312 Function Value [1] 815.6329 Gradient: [1] 9.726484e-25 -2.590048e-08 Successive iterates within tolerance. Current iterate is probably solution. iteration = 0 Step: [1] 0 0 Parameter: [1] 0.7142093 0.8328863 Function Value [1] 1045.966 Gradient: [1] -2.701254 1.837528 iteration = 6 Parameter: [1] 1.1387817 0.8237138 Function Value [1] 1045.492 Gradient: [1] -3.297258e-08 4.558449e-08 Successive iterates within tolerance. Current iterate is probably solution. iteration = 0 Step: [1] 0 0 Parameter: [1] 1.2190786 0.8233033 Function Value [1] 1044.609 Gradient: [1] 2.589381 -2.443006 iteration = 7 Parameter: [1] 0.6880674 0.8342263 Function Value [1] 1043.658 Gradient: [1] -5.457107e-08 -1.710382e-06 Successive iterates within tolerance. Current iterate is probably solution. iteration = 0 Step: [1] 0 0 Parameter: [1] 0.6867763 0.8348800 Function Value [1] 1045.262 Gradient: [1] -3.318076 2.320207 iteration = 7 Parameter: [1] 1.2547275 0.8239144 Function Value [1] 1044.541 Gradient: [1] 1.496284e-08 -6.761539e-21 Successive iterates within tolerance. Current iterate is probably solution. iteration = 0 Step: [1] 0 0 Parameter: [1] 1.2661227 0.8244986 Function Value [1] 1044.441 Gradient: [1] 2.642391 -2.502298 iteration = 1 Parameter: [1] -1509.6490716 0.8231054 Function Value [1] 826.5788 Gradient: [1] 5.044484e-19 3.488575e+01 Maximum step size exceeded 5 consecutive times. Either the function is unbounded below, becomes asymptotic to a finite value from above in some direction, or stepmx is too small. iteration = 0 Step: [1] 0 0 Parameter: [1] 0.6851074 0.8350764 Function Value [1] 1045.189 Gradient: [1] -3.359032 2.357977 iteration = 7 Parameter: [1] 1.2648770 0.8239819 Function Value [1] 1044.448 Gradient: [1] -1.484278e-08 -4.556966e-08 Last global step failed to locate a point lower than x. Either x is an approximate local minimum of the function, the function is too non-linear for this algorithm, or steptol is too large. iteration = 0 Step: [1] 0 0 Parameter: [1] 1.2679928 0.8245302 Function Value [1] 1044.427 Gradient: [1] 2.646834 -2.508039 iteration = 1 Parameter: [1] -1511.23183 0.82312 Function Value [1] 826.1627 Gradient: [1] 1.373859e-18 3.405617e+01 Maximum step size exceeded 5 consecutive times. Either the function is unbounded below, becomes asymptotic to a finite value from above in some direction, or stepmx is too small. iteration = 0 Step: [1] 0 0 Parameter: [1] 0.6849866 0.8350863 Function Value [1] 1045.184 Gradient: [1] -3.361700 2.360219 iteration = 7 Parameter: [1] 1.2654792 0.8239847 Function Value [1] 1044.442 Gradient: [1] -1.483571e-08 -2.278475e-08 Successive iterates within tolerance. Current iterate is probably solution. iteration = 0 Step: [1] 0 0 Parameter: [1] 1.268088 0.824533 Function Value [1] 1044.426 Gradient: [1] 2.647101 -2.508261 iteration = 1 Parameter: [1] -1511.3133956 0.8231232 Function Value [1] 827.7616 Gradient: [1] 1.046018e-18 3.410108e+01 Maximum step size exceeded 5 consecutive times. Either the function is unbounded below, becomes asymptotic to a finite value from above in some direction, or stepmx is too small. iteration = 0 Step: [1] 0 0 Parameter: [1] 0.6850998 0.8350804 Function Value [1] 1045.188 Gradient: [1] -3.359382 2.358456 iteration = 7 Parameter: [1] 1.2650078 0.8239839 Function Value [1] 1044.447 Gradient: [1] -1.484124e-08 2.278477e-08 Last global step failed to locate a point lower than x. Either x is an approximate local minimum of the function, the function is too non-linear for this algorithm, or steptol is too large. iteration = 0 Step: [1] 0 0 Parameter: [1] 1.2680136 0.8245307 Function Value [1] 1044.426 Gradient: [1] 2.646877 -2.508084 iteration = 1 Parameter: [1] -1511.2495885 0.8231376 Function Value [1] 827.6179 Gradient: [1] 1.473509e-19 3.477226e+01 Maximum step size exceeded 5 consecutive times. Either the function is unbounded below, becomes asymptotic to a finite value from above in some direction, or stepmx is too small. iteration = 0 Step: [1] 0 0 Parameter: [1] 0.6851144 0.8350799 Function Value [1] 1045.188 Gradient: [1] -3.359088 2.358286 iteration = 6 Parameter: [1] 1.2649473 0.8239835 Function Value [1] 1044.448 Gradient: [1] 2.968390e-08 -6.835435e-08 Successive iterates within tolerance. Current iterate is probably solution. iteration = 0 Step: [1] 0 0 Parameter: [1] 1.267994 0.824532 Function Value [1] 1044.427 Gradient: [1] 2.646853 -2.507863 iteration = 1 Parameter: [1] -1511.2336349 0.8231352 Function Value [1] 828.0258 Gradient: [1] -6.701743e-20 3.365341e+01 Maximum step size exceeded 5 consecutive times. Either the function is unbounded below, becomes asymptotic to a finite value from above in some direction, or stepmx is too small. iteration = 0 Step: [1] 0 0 Parameter: [1] 0.6851079 0.8350803 Function Value [1] 1045.188 Gradient: [1] -3.359215 2.358392 iteration = 7 Parameter: [1] 1.2649721 0.8239836 Function Value [1] 1044.447 Gradient: [1] 1.484166e-08 4.556957e-08 Last global step failed to locate a point lower than x. Either x is an approximate local minimum of the function, the function is too non-linear for this algorithm, or steptol is too large. iteration = 0 Step: [1] 0 0 Parameter: [1] 1.268001 0.824530 Function Value [1] 1044.427 Gradient: [1] 2.646867 -2.508106 iteration = 1 Parameter: [1] -1511.2385145 0.8231441 Function Value [1] 825.0983 Gradient: [1] -9.146131e-19 3.823514e+01 Maximum step size exceeded 5 consecutive times. Either the function is unbounded below, becomes asymptotic to a finite value from above in some direction, or stepmx is too small. iteration = 0 Step: [1] 0 0 Parameter: [1] 0.6849865 0.8350865 Function Value [1] 1045.184 Gradient: [1] -3.361722 2.360234 iteration = 6 Parameter: [1] 1.265490 0.823985 Function Value [1] 1044.442 Gradient: [1] -1.254955e-20 -6.835423e-08 Successive iterates within tolerance. Current iterate is probably solution. iteration = 0 Step: [1] 0 0 Parameter: [1] 1.2681041 0.8245325 Function Value [1] 1044.426 Gradient: [1] 2.647101 -2.508336 iteration = 1 Parameter: [1] -1511.3262872 0.8231335 Function Value [1] 828.0964 Gradient: [1] 2.680533e-19 3.401084e+01 Maximum step size exceeded 5 consecutive times. Either the function is unbounded below, becomes asymptotic to a finite value from above in some direction, or stepmx is too small. iteration = 0 Step: [1] 0 0 Parameter: [1] 0.6850990 0.8350805 Function Value [1] 1045.188 Gradient: [1] -3.359395 2.358484 iteration = 7 Parameter: [1] 1.2650090 0.8239838 Function Value [1] 1044.447 Gradient: [1] -4.452368e-08 -2.278478e-08 Last global step failed to locate a point lower than x. Either x is an approximate local minimum of the function, the function is too non-linear for this algorithm, or steptol is too large. iteration = 0 Step: [1] 0 0 Parameter: [1] 1.2680059 0.8245313 Function Value [1] 1044.426 Gradient: [1] 2.646886 -2.508011 iteration = 1 Parameter: [1] -1511.24342 0.82312 Function Value [1] 827.5457 Gradient: [1] -1.234659e-18 3.430500e+01 Maximum step size exceeded 5 consecutive times. Either the function is unbounded below, becomes asymptotic to a finite value from above in some direction, or stepmx is too small. iteration = 0 Step: [1] 0 0 Parameter: [1] 0.6851253 0.8350791 ==> 2369.reply.1 <== From bates@bates5.stat.wisc.edu Fri Dec 13 00:13:06 2002 Received: from bates5.stat.wisc.edu (mail@bates5.stat.wisc.edu [128.105.174.135]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id AAA28071 for ; Fri, 13 Dec 2002 00:13:06 +0100 (MET) Received: from bates by bates5.stat.wisc.edu with local (Exim 3.36 #1 (Debian)) id 18McV4-0002vX-00; Thu, 12 Dec 2002 17:12:14 -0600 To: jerome@hivnet.ubc.ca Cc: r-devel@stat.math.ethz.ch, R-bugs@biostat.ku.dk Subject: Re: [Rd] convergence problem with nlme() (PR#2369) References: <200212122243.XAA27894@pubhealth.ku.dk> From: Douglas Bates Date: 12 Dec 2002 17:12:14 -0600 In-Reply-To: <200212122243.XAA27894@pubhealth.ku.dk> Message-ID: <6radjade7l.fsf@bates5.stat.wisc.edu> Lines: 12 User-Agent: Gnus/5.0808 (Gnus v5.8.8) XEmacs/21.4 (Common Lisp) MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Sender: Douglas Bates It turns out that you learn more if you set verbose = TRUE in the call to nlme rather than setting control = nlmeControl(msVerbose = TRUE). The verbose = TRUE option shows you that the parameter values are bouncing back and forth between two regions of the parameter space, neither of which are close to the optimum. It happens that nlm will occasionally take very large steps at the beginning of the optimization resulting in unusual parameter values. This example is a difficult optimization problem. It may be possible to stabilize it somewhat by replacing the internal calls to nlm by calls to optim but that brings its own set of difficulties. ==> 2369.audit <== Sun Dec 15 08:53:57 2002 ripley foobar Sun Dec 15 08:53:57 2002 ripley moved from incoming to Add-ons ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 2374 * Subject: quantreg package - predict method for rq objects From: John Maindonald Date: Sat, 14 Dec 2002 20:46:34 +1100 (EST) -- no notes -- ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 2374 <== From mailnull@stat.math.ethz.ch Sat Dec 14 10:47:34 2002 Received: from hypatia.math.ethz.ch (root@hypatia.ethz.ch [129.132.58.23]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id KAA08577 for ; Sat, 14 Dec 2002 10:47:33 +0100 (MET) Received: from hypatia.math.ethz.ch (mailnull@localhost [127.0.0.1]) by hypatia.math.ethz.ch (8.12.6/8.12.6) with ESMTP id gBE9khgo024271 (version=TLSv1/SSLv3 cipher=EDH-RSA-DES-CBC3-SHA bits=168 verify=NO) for ; Sat, 14 Dec 2002 10:46:44 +0100 (MET) Received: (from mailnull@localhost) by hypatia.math.ethz.ch (8.12.6/8.12.6/Submit) id gBE9khMO024270 for R-bugs@biostat.ku.dk; Sat, 14 Dec 2002 10:46:43 +0100 (MET) Received: from franz.stat.wisc.edu (mail@www.omegahat.org [128.105.174.32]) by hypatia.math.ethz.ch (8.12.6/8.12.6) with ESMTP id gBE9kcgn024245 for ; Sat, 14 Dec 2002 10:46:40 +0100 (MET) Received: from anumail4.anu.edu.au ([150.203.2.44] helo=anu.edu.au) by franz.stat.wisc.edu with esmtp (Exim 3.35 #1 (Debian)) id 18N8sY-0008MD-00 for ; Sat, 14 Dec 2002 03:46:38 -0600 Received: from leonard.anu.edu.au (leonard.anu.edu.au [150.203.2.15]) by anu.edu.au (8.12.4/8.12.4) with ESMTP id gBE9kZBo015610 for ; Sat, 14 Dec 2002 20:46:35 +1100 (EST) Received: (from u9801539@localhost) by leonard.anu.edu.au (8.11.6/8.11.6) id gBE9kYl17197 for r-bugs@r-project.org; Sat, 14 Dec 2002 20:46:34 +1100 (EST) Date: Sat, 14 Dec 2002 20:46:34 +1100 (EST) From: John Maindonald Message-Id: <200212140946.gBE9kYl17197@leonard.anu.edu.au> To: r-bugs@r-project.org Subject: quantreg package - predict method for rq objects X-Sender: u9801539@leonard.anu.edu.au X-Sender-Domain: leonard.anu.edu.au X-Spam-Score: (0.4) X-Spam-Tests: FROM_ENDS_IN_NUMS X-Scanned-By: MIMEDefang 2.15 (www dot roaringpenguin dot com slash mimedefang) X-Virus-Scanned: by amavisd-milter (http://amavis.org/) X-Virus-Scanned: by amavisd-milter (http://amavis.org/) X-Spam-Status: No, hits=1.7 required=5.0 tests=FROM_ENDS_IN_NUMS,SPAM_PHRASE_00_01 version=2.43 X-Spam-Level: * In help(rq.object), under Methods, it says: " The `"rq"' class of objects has methods for the following generic functions: `coef', `effects' , `formula' , `labels' , `model.frame' , `model.matrix' , `plot' , `predict' , `print' , `print.summary' , `residuals' , `summary' " There seems to be no predict method: > y <- rnorm(500) > u <- rq(y ~ 1, tau=.75) Warning message: Solution may be nonunique in: rq.fit.br(x, y, tau = tau, ...) > predict(u) Error in predict(u) : no applicable method for "predict" [I have looked both at the Windows and at the Darwin/X11 versions.] John Maindonald email : john.maindonald@anu.edu.au Centre for Bioinformation Science, phone : (6125)3473 c/o MSI, fax : (6125)5549 John Dedman Mathematical Sciences Building (Building 27) Australian National University Canberra ACT 0200 Australia ==> 2374.audit <== Sun Dec 15 08:53:15 2002 ripley moved from incoming to Add-ons ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 2386 * Subject: ace and avas in acepack From: Frank E Harrell Jr Date: Fri, 20 Dec 2002 12:59:07 -0500 -- no notes -- ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 2386 <== From fharrell@virginia.edu Fri Dec 20 18:59:26 2002 Received: from cgatepro-2.mail.virginia.edu (neon.mail.Virginia.EDU [128.143.2.220]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id SAA06856 for ; Fri, 20 Dec 2002 18:59:26 +0100 (MET) Received: from [128.143.108.90] (HELO fharrell.biostat) by cgatepro-2.mail.virginia.edu (CommuniGate Pro SMTP 3.5.9) with SMTP id 29368136; Fri, 20 Dec 2002 12:58:38 -0500 Date: Fri, 20 Dec 2002 12:59:07 -0500 From: Frank E Harrell Jr To: r-bugs@biostat.ku.dk Cc: Thomas Lumley Subject: ace and avas in acepack Message-Id: <20021220125907.57c1ad10.fharrell@virginia.edu> Organization: University of Virginia X-Mailer: Sylpheed version 0.8.2 (GTK+ 1.2.10; i386-redhat-linux-gnu) Mime-Version: 1.0 Content-Type: text/plain; charset=US-ASCII Content-Transfer-Encoding: 7bit Under platform i686-pc-linux-gnu arch i686 os linux-gnu system i686, linux-gnu status major 1 minor 6.1 year 2002 month 11 day 01 language R and using acepack 1.3-2, ace() and avas() have checks such as if(circ[i] < 0 || circ[i] > nrow(x)) { if(mon[i] < 0 || mon[i] > nrow(x)) { if(lin[i] < 0 || lin[i] > nrow(x)) { etc. I believe that nrow should be ncol. Also, in ace() there is the follwing check: if (cat[i] == 0) { cat("response spec can only be lin or ordered (default)") return() } whereas I believe that the ace algorithm does allow for categorical response variables. -- Frank E Harrell Jr Prof. of Biostatistics & Statistics Div. of Biostatistics & Epidem. Dept. of Health Evaluation Sciences U. Virginia School of Medicine http://hesweb1.med.virginia.edu/biostat ==> 2386.audit <== Fri Dec 27 09:19:24 2002 ripley moved from incoming to Add-ons ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 2480 * Subject: bug in CrossTable (package:gregmisc) From: John_Hendrickx@yahoo.com Date: Tue, 21 Jan 2003 13:51:01 +0100 (MET) -- no notes -- ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 2480 <== From John_Hendrickx@yahoo.com Tue Jan 21 13:51:02 2003 Received: from blueberry (blueberry [192.38.19.4]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id NAA17686 for ; Tue, 21 Jan 2003 13:51:01 +0100 (MET) Date: Tue, 21 Jan 2003 13:51:01 +0100 (MET) From: John_Hendrickx@yahoo.com Message-Id: <200301211251.NAA17686@pubhealth.ku.dk> To: R-bugs@biostat.ku.dk Subject: bug in CrossTable (package:gregmisc) Full_Name: John Hendrickx Version: 1.6.0 OS: Windows 98 Submission from: (NULL) (137.224.174.216) CrossTable in the "gregmisc" package fails when the fisher.exact test produces an error (I suspect this is because the number of cases is too large). This can be fixed using "FTt <- try(fisher.test(t, alternative = "two.sided"))" or by making the test optional. bugtab <- array(c( 24, 31, 2, 50, 4, 101, 16, 130, 0, 29, 3, 52, 3, 39, 11, 210)) dim(bugtab)<-c(4,4) bugtab<-t(bugtab) rownames(bugtab)<-c("a","b","c","d") colnames(bugtab)<-c("A","B","C","D") bugtab thisworks<-trunc(bugtab/10) thisworks library(gregmisc) # fisher.exact test works fine on small tables CrossTable(thisworks) # I suspect the bugtab table is too large CrossTable(bugtab) ==> 2480.audit <== Tue Jan 21 14:39:22 2003 ripley moved from incoming to Add-ons ==> 2480.followup.1 <== From mschwartz@medanalytics.com Wed Jan 22 05:23:39 2003 Received: from mail12.atl.registeredsite.com (nobody@mail12.atl.registeredsite.com [64.224.219.86]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id FAA22283 for ; Wed, 22 Jan 2003 05:23:38 +0100 (MET) Received: from mail.medanalytics.com (mail.medanalytics.com [64.225.164.192]) by mail12.atl.registeredsite.com (8.12.2/8.12.6) with ESMTP id h0M4Mn0s024550; Tue, 21 Jan 2003 23:22:49 -0500 Received: from MARC [64.225.164.192] by mail.medanalytics.com with ESMTP (SMTPD32-6.06) id AC96524D00F2; Tue, 21 Jan 2003 23:22:46 -0500 Reply-To: From: "Marc Schwartz" To: , Cc: Subject: RE: [Rd] bug in CrossTable (package:gregmisc) (PR#2480) Date: Tue, 21 Jan 2003 22:22:45 -0600 Organization: MedAnalytics, Inc. Message-ID: <009801c2c1cd$eaebe5e0$0201a8c0@MARC> MIME-Version: 1.0 Content-Type: text/plain; charset="US-ASCII" X-Priority: 3 (Normal) X-MSMail-Priority: Normal X-Mailer: Microsoft Outlook, Build 10.0.4510 Importance: Normal In-Reply-To: <200301211251.NAA17691@pubhealth.ku.dk> X-MimeOLE: Produced By Microsoft MimeOLE V6.00.2800.1106 Content-Transfer-Encoding: 8bit X-MIME-Autoconverted: from quoted-printable to 8bit by pubhealth.ku.dk id FAA22283 > -----Original Message----- > From: r-devel-admin@stat.math.ethz.ch > [mailto:r-devel-admin@stat.math.ethz.ch] On Behalf Of > John_Hendrickx@yahoo.com > Sent: Tuesday, January 21, 2003 6:51 AM > To: r-devel@stat.math.ethz.ch > Cc: R-bugs@biostat.ku.dk > Subject: [Rd] bug in CrossTable (package:gregmisc) (PR#2480) > > > Full_Name: John Hendrickx > Version: 1.6.0 > OS: Windows 98 > Submission from: (NULL) (137.224.174.216) > > > CrossTable in the "gregmisc" package fails when the > fisher.exact test produces an error (I suspect this is > because the number of cases is too large). This can be fixed > using "FTt <- try(fisher.test(t, alternative = "two.sided"))" > or by making the test optional. > > bugtab <- array(c( > 24, 31, 2, 50, > 4, 101, 16, 130, > 0, 29, 3, 52, > 3, 39, 11, 210)) > dim(bugtab)<-c(4,4) > bugtab<-t(bugtab) > rownames(bugtab)<-c("a","b","c","d") > colnames(bugtab)<-c("A","B","C","D") > bugtab > > thisworks<-trunc(bugtab/10) > thisworks > > library(gregmisc) > # fisher.exact test works fine on small tables > CrossTable(thisworks) > # I suspect the bugtab table is too large > CrossTable(bugtab) John, Thank you for bringing this to my attention. I was out of town all day on a business trip and will look at this issue further with fresh eyes in the morning. I'll get an update to Greg as soon as I have a chance to incorporate a fix. I also have some pending updates for barplot2(), so I will coordinate the updates with Greg to conserve his time. Best regards, Marc Schwartz ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ * PR# 2499 * Subject: survival bug? From: Pius Korner Date: Mon, 27 Jan 2003 17:28:11 +0100 --Needs to resubmit in more useful format. ~~~~~~~~~~ original reports follow ~~~~~~~~~~ ==> 2499 <== From korner@eco.umnw.ethz.ch Mon Jan 27 17:28:27 2003 Received: from xsmtp.ethz.ch (xsmtp.ethz.ch [129.132.97.6]) by pubhealth.ku.dk (8.9.3/8.9.1) with ESMTP id RAA07671 for ; Mon, 27 Jan 2003 17:28:26 +0100 (MET) Received: from xfe2.d.ethz.ch ([192.168.36.11]) by xsmtp.ethz.ch with Microsoft SMTPSVC(5.0.2195.5329); Mon, 27 Jan 2003 17:27:39 +0100 Received: from eco.umnw.ethz.ch ([129.132.15.199]) by xfe2.d.ethz.ch with Microsoft SMTPSVC(5.0.2195.5329); Mon, 27 Jan 2003 17:27:38 +0100 Date: Mon, 27 Jan 2003 17:28:11 +0100 Mime-Version: 1.0 (Apple Message framework v548) Content-Type: multipart/mixed; boundary=Apple-Mail-27-953181986 Subject: survival bug? From: Pius Korner To: r-bugs@biostat.ku.dk Message-Id: <5397D74A-3214-11D7-912D-000A27D7DD96@eco.umnw.ethz.ch> X-Mailer: Apple Mail (2.548) X-OriginalArrivalTime: 27 Jan 2003 16:27:38.0955 (UTC) FILETIME=[01FFB1B0:01C2C621] --Apple-Mail-27-953181986 Content-Transfer-Encoding: 7bit Content-Type: text/plain; charset=US-ASCII; format=flowed a possible bug with survival analysis - either in R or in SPSS... find more details in bug.doc, and the data in bug.txt best Pius Korner --Apple-Mail-27-953181986 Content-Disposition: attachment Content-Type: multipart/appledouble; boundary=Apple-Mail-28-953181987 --Apple-Mail-28-953181987 Content-Disposition: attachment; filename=bug.doc Content-Transfer-Encoding: base64 Content-Type: application/applefile; name="bug.doc" AAUWBwACAAAAAAAAAAAAAAAAAAAAAAAAAAMAAAAJAAAAPgAAAAoAAAADAAAASAAAAAcAAAACAAAA TwAAAR5XOEJOTVNXRAAAYnVnLmRvYwAAAQAAAAEAAAAAAAAAAB4AAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAEAAAAB AAAAAAAAAAAeAAAAAAAAAAAAHAAe//8= --Apple-Mail-28-953181987 Content-Disposition: attachment; filename=bug.doc Content-Transfer-Encoding: base64 Content-Type: application/msword; x-mac-creator=4D535744; x-unix-mode=0644; x-mac-type=5738424E; name="bug.doc" 0M8R4KGxGuEAAAAAAAAAAAAAAAAAAAAAPgADAP7/CQAGAAAAAAAAAAAAAAABAAAAJQAAAAAAAAAA EAAAJwAAAAEAAAD+////AAAAACQAAAD///////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////s pcEAsUAJBAAA8BK/AAAAAAABEQABAAEABgAAUwkAAA4AamJqYpgKmAoAAAAAAAAAAAAAAAAAAAAA AAAJBBYAHhYAAPJoAQDyaAEAUwMAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAD//w8AAAAA AAAAAAD//w8AAAAAAAAAAAD//w8AAAAAAAAAAAAAAAAAAAAAAGwAAAAAAKIAAAAAAAAAogAAAKIA AAAAAAAAogAAAAAAAACiAAAAAAAAAKIAAAAAAAAAogAAABQAAAAAAAAAAAAAANIAAAAAAAAA0gAA AAAAAADSAAAAAAAAANIAAAAAAAAA0gAAAAwAAADeAAAAFAAAANIAAAAAAAAAswEAAGACAAD+AAAA AAAAAP4AAAAAAAAA/gAAAAAAAAD+AAAAAAAAAP4AAAAAAAAA/gAAAAAAAAD+AAAAAAAAAP4AAAAA AAAAcAEAAAIAAAByAQAAAAAAAHIBAAAAAAAAcgEAAAAAAAByAQAAAAAAAHIBAAAAAAAAcgEAACwA AAATBAAAIAIAADMGAABiAAAAngEAABUAAAAAAAAAAAAAAAAAAAAAAAAAogAAAAAAAAD+AAAAAAAA AAAAAAAAAAAAAAAAAAAAAAD+AAAAAAAAAP4AAAAAAAAA/gAAAAAAAAD+AAAAAAAAAJ4BAAAAAAAA /gAAAAAAAACiAAAAAAAAAKIAAAAAAAAA/gAAAAAAAAAAAAAAAAAAAP4AAAAAAAAA/gAAAAAAAAD+ AAAAAAAAAP4AAAAAAAAA/gAAAAAAAAD+AAAAAAAAAKIAAAAAAAAA/gAAAAAAAACiAAAAAAAAAP4A AAAAAAAAcAEAAAAAAAAAAAAAAAAAAP4AAAAAAAAAtgAAAA4AAADEAAAADgAAAKIAAAAAAAAAogAA AAAAAACiAAAAAAAAAKIAAAAAAAAA/gAAAAAAAABwAQAAAAAAAP4AAAByAAAA/gAAAAAAAAAAAAAA AAAAAHABAAAAAAAAogAAAAAAAACiAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAcAEAAAAAAAAAAAAAAAAAAPIAAAAMAAAAiw1bugAA AADSAAAAAAAAANIAAAAAAAAA/gAAAAAAAABwAQAAAAAAAAAAAAAAAAAAcAEAAAAAAACzAQAAAAAA ALMBAAAAAAAAcAEAAAAAAACVBgAAAAAAAP4AAAAAAAAAlQYAAAAAAABwAQAAAAAAAP4AAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABADDAA8A5QAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAABEZWFy IFINDWlzIHRoYXQgYSBidWc/PzoNDVItY29kZToNIyBkYXRhIHByZXBhcmF0aW9uDWQuZGF0YTwt cmVhZC50YWJsZSgiYnVnLnR4dCIsaGVhZGVyPVQpDWQuZGF0YSRmYWN0MTwtYXMuZmFjdG9yKGQu RjEkZmFjdDEpDWQuZGF0YSRmYWN0MjwtYXMuZmFjdG9yKGQuRjEkZmFjdDIpDWxpYnJhcnkoc3Vy dml2YWwpDSMgQ294IFByb3BvcnRpb25hbCBIYXphcmQNci5jb3g8LWNveHBoKFN1cnYodGltZSxk ZWx0YSl+ZmFjdDErZmFjdDIrY292YXIsZGF0YT1kLmRhdGEpDWRyb3AxKHIuY294LHRlc3Q9IkNo aXNxIikNDW91dHB1dDoNU2luZ2xlIHRlcm0gZGVsZXRpb25zDQ1Nb2RlbDoNU3Vydih0aW1lLCBk ZWx0YSkgfiBmYWN0MSArIGZhY3QyICsgY292YXINICAgICAgIERmICAgIEFJQyAgICBMUlQgICBQ cihDaGkpICAgIA08bm9uZT4gICAgODI3LjkwICAgICAgICAgICAgICAgICAgICAgDWZhY3QxICAg NyA4MjguNDkgIDE0LjU4ICAgMC4wNDE3MiAqICANZmFjdDIgICA0IDg0My41NCAgMjMuNjMgOS40 NTVlLTA1ICoqKg1jb3ZhciAgIDEgODI1Ljk4ICAgMC4wOCAgIDAuNzc3NjkgICAgDS0tLQ1TaWdu aWYuIGNvZGVzOiAgMCBgKioqJyAwLjAwMSBgKionIDAuMDEgYConIDAuMDUgYC4nIDAuMSBgICcg MQ0NDWlmIEkgZG8gdGhlIHNhbWUgb24gc3BzcywgdGhlIHAtdmFsdWUgZm9yIGZhY3QxIGlzIDAu MDIsIGZvciBmYWN0MiwgaXRzIDAuMDAxLiBJcyB0aGF0IGEgYnVnIG9yIGp1c3QgZGlmZmVyZW50 IHByb2NlZHVyZXM/DQ1hbHNvLCBpZiBOQSdzIGFyZSBpbiB0aGUgZmlsZSwgdGhlIHAtdmFsdWVz IHNlZW0gdG8gYmUgY29tcGxldGVseSBpbXBvc3NpYmxlLi4uDQ1iZXN0DVBpdXMNDQAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGAAAiBgAA JAYAADQGAABQBgAAtgYAAM8GAAAnBwAAMAcAAFIJAABTCQAAAPnx+eri6gDZ1AAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAIT0oHAFFKBwAAEUIqAU9KBwBRSgcAcGgAAAAADzUIgU9KBgBQSgAA UUoGAAxPSgUAUEoAAFFKBQAADzUIgU9KBABQSgAAUUoEAAxPSgMAUEoAAFFKAwAKAAYAAAcGAAAI BgAAGQYAABoGAAAiBgAANQYAAFwGAACABgAApAYAALYGAADQBgAADQcAACcHAAAoBwAAMAcAAEYH AABHBwAATgcAAHgHAACeBwAAxAcAAOoHAAAQCAAANggAADoIAAB5CAAA/QAAAAAAAAAAAAAAAP0A AAAAAAAAAAAAAAD9AAAAAAAAAAAAAAAA/QAAAAAAAAAAAAAAAP0AAAAAAAAAAAAAAADhAAAAAAAA AAAAAAAA4QAAAAAAAAAAAAAAAOEAAAAAAAAAAAAAAADhAAAAAAAAAAAAAAAA4QAAAAAAAAAAAAAA AOEAAAAAAAAAAAAAAADhAAAAAAAAAAAAAAAA4QAAAAAAAAAAAAAAAP0AAAAAAAAAAAAAAAD9AAAA AAAAAAAAAAAA/QAAAAAAAAAAAAAAAP0AAAAAAAAAAAAAAAD9AAAAAAAAAAAAAAAA/QAAAAAAAAAA AAAAAP0AAAAAAAAAAAAAAAD9AAAAAAAAAAAAAAAA/QAAAAAAAAAAAAAAAP0AAAAAAAAAAAAAAAD9 AAAAAAAAAAAAAAAA/QAAAAAAAAAAAAAAAP0AAAAAAAAAAAAAAAAAAAAAHAAAMSQANyQAOCQASCQA DcYmAAzQAqAFcAhACxAO4BCwE4AWUBkgHPAewCEAAAAAAAAAAAAAAAAAAQAAABoABgAABwYAAAgG AAAZBgAAGgYAACIGAAA1BgAAXAYAAIAGAACkBgAAtgYAANAGAAANBwAAJwcAACgHAAAwBwAARgcA AEcHAABOBwAAeAcAAJ4HAADEBwAA6gcAABAIAAA2CAAAOggAAHkIAAB6CAAAewgAAPYIAAD3CAAA RwkAAEgJAABNCQAAUgkAAFMJAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAI3kIAAB6CAAAewgA APYIAAD3CAAARwkAAEgJAABNCQAAUgkAAFMJAAD9AAAAAAAAAAAAAAAA/QAAAAAAAAAAAAAAAP0A AAAAAAAAAAAAAAD9AAAAAAAAAAAAAAAA/QAAAAAAAAAAAAAAAP0AAAAAAAAAAAAAAAD9AAAAAAAA AAAAAAAA/QAAAAAAAAAAAAAAAP0AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAEAAAAJHAAfsNAvILDgPSGw CAcisAgHI5CgBSSQoAUlsAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAASAA8ACgABAGkADwAC AAgACAAIADQAAEDx/wIANAAAAAYATgBvAHIAbQBhAGwAAAACAAAAFABDShgAT0oIAFBKCABRSggA bUgJBAAAAAAAAAAAAAAAAAAAAAAAADwAQUDy/6EAPAAAABYARABlAGYAYQB1AGwAdAAgAFAAYQBy AGEAZwByAGEAcABoACAARgBvAG4AdAAAAAAAAAAAAAAAAAAAAAAAUwMAAAUAABYAAAAA/////wEA AAAEIP//AQCgepkAAAAAAFMDAAAAAAAAAAAABgAAUwkAAAcAAAAABgAAeQgAAFMJAAAIAAAACgAA AAAGAABTCQAACQAAAP//AgAAAAsAcABpAHUAcwAgAGsAbwByAG4AZQByACkAVABvAHAAbwBsAGkA bgBvACAASABEACAAOgBVAHMAZQByAHMAOgBrAG8AcgBuAGUAcgA6AEQAZQBzAGsAdABvAHAAOgBi AHUAZwAuAGQAbwBjAP9AA4ABAFIDAABSAwAA3GOyAwEAAQBSAwAAAAAAAHgCAAAAAAAAAwADACMD owICEAAAAAAAAABTAwAAUAAADABAAAAJAAAARwaQAQAAAAICBgMFBAUCAwAAAAMAAAAAAAAAAAAA AAABAAAAAAAAAFQAaQBtAGUAcwAgAE4AZQB3ACAAUgBvAG0AYQBuAAAANQaQAQIAAAIABQAAAAAA AAAAAAAQAAAAAAAAAAAAAAAAAACAAAAAAFMAeQBtAGIAbwBsAAAAMwaQAQAAAAILBgQCAgICAgAA AAMAAAAAAAAAAAAAAAABAAAAAAAAAEEAcgBpAGEAbAAAADsGkAEAAAACAAUAAAAAAAAAAAADAAAA AAAAAAAAAAAAAQAAAAAAAABIAGUAbAB2AGUAdABpAGMAYQAAAFEgkAFNDwAAAAAAAAAAAAAAAAAD AAAAAAAAAAAAAAAAAQAAAAAAAABIAGUAbAB2AGUAdABpAGMAYQAtAEIAbwBsAGQAAABUAGkAbQBl AHMAAABNAJABTQ0AAAAAAAAAAAAAAAAAAwAAAAAAAAAAAAAAAAEAAAAAAAAATAB1AGMAaQBkAGEA RwByAGEAbgBkAGUAAABUAGkAbQBlAHMAAABXAJABTRIAAAAAAAAAAAAAAAAAAwAAAAAAAAAAAAAA AAEAAAAAAAAATAB1AGMAaQBkAGEARwByAGEAbgBkAGUALQBCAG8AbABkAAAAVABpAG0AZQBzAAAA QQaQAQAAAAILBQYCAgIDAgAAAAMAAAAAAAAAAAAAAAABAAAAAAAAAEEAcgBpAGEAbAAgAE4AYQBy AHIAbwB3AAAAMwaQAQAAAAIABQAAAAAAAAAAAAMAAAAAAAAAAAAAAAABAAAAAAAAAFQAaQBtAGUA cwAAACAABAAxCIwYAPDQAgAAaAEAAAAAVNxxJlfccSYAAAAAAQADAAAAAAAAAAAAAAABAAEAAAAE AAMQAQAAAAAAAAAAAAAAAQABAAAAAQAAAAAAAADpAwDwEAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAClBsAHtAC0AIAAEjAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAgAAAHICAAAAAAgAAEAA 8BAACAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAD//xIAAAAAAAAABgBEAGUAYQByACAA UgAAAAAAAAALAHAAaQB1AHMAIABrAG8AcgBuAGUAcgALAHAAaQB1AHMAIABrAG8AcgBuAGUAcgAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA AAAAAAAAAAAAAAAAAA