Re: [Rd] Calculation of e^{z^2/2} for a normal deviate z
On 22/06/2019 00:58, jing hua zhao wrote:
Hi Peter, Rui, Chrstophe and Gabriel,
Thanks for your inputs -- the use of qnorm(., log=TRUE) is a good point
Another approach could be simply to note that a function defined as
f(p)=exp(-z(p)^2/2) is regular around p=0 with f(0)=0.
It has roughly the shape of p*(2-p) for p \in [0; 1]. So we can
calculate let say f(10^-10) with sufficient precision using Rmpfr and
then use a linear approximation for p from [0, 10^-10]. After that a
simple inverse gives us e^(z*z/2).
Serguei.
in line with pnorm with which we devised log(p) as
log(2) + pnorm(-abs(z), lower.tail = TRUE, log.p = TRUE)
that could do really really well for large z compared to Rmpfr. Maybe I am
asking too much since
z <-2
Rmpfr::format(2*pnorm(mpfr(-abs(z),100),lower.tail=TRUE,log.p=FALSE))
[1] "1.660579603192917090365313727164e-86858901"
already gives a rarely seen small p value. I gather I also need a multiple
precision exp() and their sum since exp(z^2/2) is also a Bayes Factor so I get
log(x_i )/sum_i log(x_i) instead. To this point, I am obliged to clarify, see
https://statgen.github.io/gwas-credible-sets/method/locuszoom-credible-sets.pdf.
I agree many feel geneticists go to far with small p values which I would have
difficulty to argue againston the other hand it is also expected to see these
in a non-genetic context. For instance the Framingham study was established in
1948 just got $34m for six years on phenotypewide association which we would be
interesting to see.
Best wishes,
Jing Hua
From: peter dalgaard
Sent: 21 June 2019 16:24
To: jing hua zhao
Cc: Rui Barradas; r-devel@r-project.org
Subject: Re: [Rd] Calculation of e^{z^2/2} for a normal deviate z
You may want to look into using the log option to qnorm
e.g., in round figures:
log(1e-300)
[1] -690.7755
qnorm(-691, log=TRUE)
[1] -37.05315
exp(37^2/2)
[1] 1.881797e+297
exp(-37^2/2)
[1] 5.314068e-298
Notice that floating point representation cuts out at 1e+/-308 or so. If you
want to go outside that range, you may need explicit manipulation of the log
values. qnorm() itself seems quite happy with much smaller values:
qnorm(-5000, log=TRUE)
[1] -99.94475
-pd
On 21 Jun 2019, at 17:11 , jing hua zhao wrote:
Dear Rui,
Thanks for your quick reply -- this allows me to see the bottom of this. I was
hoping we could have a handle of those p in genmoics such as 1e-300 or smaller.
Best wishes,
Jing Hua
From: Rui Barradas
Sent: 21 June 2019 15:03
To: jing hua zhao; r-devel@r-project.org
Subject: Re: [Rd] Calculation of e^{z^2/2} for a normal deviate z
Hello,
Well, try it:
p <- .Machine$double.eps^seq(0.5, 1, by = 0.05)
z <- qnorm(p/2)
pnorm(z)
# [1] 7.450581e-09 1.22e-09 2.026908e-10 3.343152e-11 5.514145e-12
# [6] 9.094947e-13 1.500107e-13 2.474254e-14 4.080996e-15 6.731134e-16
#[11] 1.110223e-16
p/2
# [1] 7.450581e-09 1.22e-09 2.026908e-10 3.343152e-11 5.514145e-12
# [6] 9.094947e-13 1.500107e-13 2.474254e-14 4.080996e-15 6.731134e-16
#[11] 1.110223e-16
exp(z*z/2)
# [1] 9.184907e+06 5.301421e+07 3.073154e+08 1.787931e+09 1.043417e+10
# [6] 6.105491e+10 3.580873e+11 2.104460e+12 1.239008e+13 7.306423e+13
#[11] 4.314798e+14
p is the smallest possible such that 1 + p != 1 and I couldn't find
anything to worry about.
R version 3.6.0 (2019-04-26)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 19.04
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.8.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.8.0
locale:
[1] LC_CTYPE=pt_PT.UTF-8 LC_NUMERIC=C
[3] LC_TIME=pt_PT.UTF-8LC_COLLATE=pt_PT.UTF-8
[5] LC_MONETARY=pt_PT.UTF-8LC_MESSAGES=pt_PT.UTF-8
[7] LC_PAPER=pt_PT.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=pt_PT.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods
[7] base
other attached packages:
[many packages loaded]
Hope this helps,
Rui Barradas
�s 15:24 de 21/06/19, jing hua zhao escreveu:
Dear R-developers,
I am keen to calculate exp(z*z/2) with z=qnorm(p/2) and p is very small. I
wonder if anyone has experience with this?
Thanks very much in advance,
Jing Hua
[[alternative HTML version deleted]]
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
[[alternative HTML version deleted]]
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd@cbs.dk Priv: pda...@gmail.com
[[alternative HTML version deleted]]
__
R-de
Re: [Rd] Calculation of e^{z^2/2} for a normal deviate z
> William Dunlap via R-devel
> on Sun, 23 Jun 2019 10:34:47 -0700 writes:
> William Dunlap via R-devel
> on Sun, 23 Jun 2019 10:34:47 -0700 writes:
> include/Rmath.h declares a set of 'logspace' functions for use at the C
> level. I don't think there are core R functions that call them.
> /* Compute the log of a sum or difference from logs of terms, i.e.,
> *
> * log (exp (logx) + exp (logy))
> * or log (exp (logx) - exp (logy))
> *
> * without causing overflows or throwing away too much accuracy:
> */
> double Rf_logspace_add(double logx, double logy);
> double Rf_logspace_sub(double logx, double logy);
> double Rf_logspace_sum(const double *logx, int nx);
> Bill Dunlap
> TIBCO Software
> wdunlap tibco.com
Yes, indeed, thank you, Bill!
But they *have* been in use by core R functions for a long time
in pgamma, pbeta and related functions.
[and I have had changes in *hyper.c where logspace_add() is used too,
for several years now (since 2015) but I no longer know which
concrete accuracy problem that addresses, so have not yet committed it]
Martin Maechler
ETH Zurich and R Core Team
> On Sun, Jun 23, 2019 at 1:40 AM Ben Bolker wrote:
>>
>> I agree with many the sentiments about the wisdom of computing very
>> small p-values (although the example below may win some kind of a prize:
>> I've seen people talking about p-values of the order of 10^(-2000), but
>> never 10^(-(10^8)) !). That said, there are a several tricks for
>> getting more reasonable sums of very small probabilities. The first is
>> to scale the p-values by dividing the *largest* of the probabilities,
>> then do the (p/sum(p)) computation, then multiply the result (I'm sure
>> this is described/documented somewhere). More generally, there are
>> methods for computing sums on the log scale, e.g.
>>
>>
>>
https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.misc.logsumexp.html
>>
>> I don't know where this has been implemented in the R ecosystem, but
>> this sort of computation is the basis of the "Brobdingnag" package for
>> operating on very large ("Brobdingnagian") and very small
>> ("Lilliputian") numbers.
>>
>>
>> On 2019-06-21 6:58 p.m., jing hua zhao wrote:
>> > Hi Peter, Rui, Chrstophe and Gabriel,
>> >
>> > Thanks for your inputs -- the use of qnorm(., log=TRUE) is a good
point
>> in line with pnorm with which we devised log(p) as
>> >
>> > log(2) + pnorm(-abs(z), lower.tail = TRUE, log.p = TRUE)
>> >
>> > that could do really really well for large z compared to Rmpfr. Maybe I
>> am asking too much since
>> >
>> > z <-2
>> >> Rmpfr::format(2*pnorm(mpfr(-abs(z),100),lower.tail=TRUE,log.p=FALSE))
>> > [1] "1.660579603192917090365313727164e-86858901"
>> >
>> > already gives a rarely seen small p value. I gather I also need a
>> multiple precision exp() and their sum since exp(z^2/2) is also a Bayes
>> Factor so I get log(x_i )/sum_i log(x_i) instead. To this point, I am
>> obliged to clarify, see
>>
https://statgen.github.io/gwas-credible-sets/method/locuszoom-credible-sets.pdf
>> .
>> >
>> > I agree many feel geneticists go to far with small p values which I
>> would have difficulty to argue againston the other hand it is also
expected
>> to see these in a non-genetic context. For instance the Framingham study
>> was established in 1948 just got $34m for six years on phenotypewide
>> association which we would be interesting to see.
>> >
>> > Best wishes,
>> >
>> >
>> > Jing Hua
>> >
>> >
>> >
>> > From: peter dalgaard
>> > Sent: 21 June 2019 16:24
>> > To: jing hua zhao
>> > Cc: Rui Barradas; r-devel@r-project.org
>> > Subject: Re: [Rd] Calculation of e^{z^2/2} for a normal deviate z
>> >
>> > You may want to look into using the log option to qnorm
>> >
>> > e.g., in round figures:
>> >
>> >> log(1e-300)
>> > [1] -690.7755
>> >> qnorm(-691, log=TRUE)
>> > [1] -37.05315
>> >> exp(37^2/2)
>> > [1] 1.881797e+297
>> >> exp(-37^2/2)
>> > [1] 5.314068e-298
>> >
>> > Notice that floating point representation cuts out at 1e+/-308 or so.
If
>> you want to go outside that range, you may need explicit manipulation of
>> the log values. qnorm() itself seems quite happy with much smaller
values:
>> >
>> >> qnorm(-5000, log=TRUE)
>> > [1] -99.94475
>> >
>> > -pd
>> >
>> >> On 21 Jun 2019, at 17:11 , jing hua zhao
>> wrote:
>> >>
>> >> Dear Rui,
>> >>
>> >> Thanks for your quick reply -- this allows me to see the bottom of
>> this. I was hoping we could have a handle of those p in genmoics such as
>> 1e-300 or smal
Re: [Rd] Calculation of e^{z^2/2} for a normal deviate z
Hi All,
Thanks for all your comments which allows me to appreciate more of these in
Python and R.
I just came across the matrixStats package,
## EXAMPLE #1
lx <- c(1000.01, 1000.02)
y0 <- log(sum(exp(lx)))
print(y0) ## Inf
y1 <- logSumExp(lx)
print(y1) ## 1000.708
and
> ly <- lx*10
> ly
[1] 11000 12000
> y1 <- logSumExp(ly)
> print(y1)
[1] 12000
> logSumExp
function (lx, idxs = NULL, na.rm = FALSE, ...)
{
has_na <- TRUE
.Call(C_logSumExp, as.numeric(lx), idxs, as.logical(na.rm),
has_na)
}
Maybe this is rather close?
Best wishes,
Jing Hua
From: R-devel on behalf of Martin Maechler
Sent: 24 June 2019 08:37
To: William Dunlap
Cc: r-devel@r-project.org
Subject: Re: [Rd] Calculation of e^{z^2/2} for a normal deviate z
> William Dunlap via R-devel
> on Sun, 23 Jun 2019 10:34:47 -0700 writes:
> William Dunlap via R-devel
> on Sun, 23 Jun 2019 10:34:47 -0700 writes:
> include/Rmath.h declares a set of 'logspace' functions for use at the C
> level. I don't think there are core R functions that call them.
> /* Compute the log of a sum or difference from logs of terms, i.e.,
> *
> * log (exp (logx) + exp (logy))
> * or log (exp (logx) - exp (logy))
> *
> * without causing overflows or throwing away too much accuracy:
> */
> double Rf_logspace_add(double logx, double logy);
> double Rf_logspace_sub(double logx, double logy);
> double Rf_logspace_sum(const double *logx, int nx);
> Bill Dunlap
> TIBCO Software
> wdunlap tibco.com
Yes, indeed, thank you, Bill!
But they *have* been in use by core R functions for a long time
in pgamma, pbeta and related functions.
[and I have had changes in *hyper.c where logspace_add() is used too,
for several years now (since 2015) but I no longer know which
concrete accuracy problem that addresses, so have not yet committed it]
Martin Maechler
ETH Zurich and R Core Team
> On Sun, Jun 23, 2019 at 1:40 AM Ben Bolker wrote:
>>
>> I agree with many the sentiments about the wisdom of computing very
>> small p-values (although the example below may win some kind of a prize:
>> I've seen people talking about p-values of the order of 10^(-2000), but
>> never 10^(-(10^8)) !). That said, there are a several tricks for
>> getting more reasonable sums of very small probabilities. The first is
>> to scale the p-values by dividing the *largest* of the probabilities,
>> then do the (p/sum(p)) computation, then multiply the result (I'm sure
>> this is described/documented somewhere). More generally, there are
>> methods for computing sums on the log scale, e.g.
>>
>>
>>
https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.misc.logsumexp.html
>>
>> I don't know where this has been implemented in the R ecosystem, but
>> this sort of computation is the basis of the "Brobdingnag" package for
>> operating on very large ("Brobdingnagian") and very small
>> ("Lilliputian") numbers.
>>
>>
>> On 2019-06-21 6:58 p.m., jing hua zhao wrote:
>> > Hi Peter, Rui, Chrstophe and Gabriel,
>> >
>> > Thanks for your inputs -- the use of qnorm(., log=TRUE) is a good
point
>> in line with pnorm with which we devised log(p) as
>> >
>> > log(2) + pnorm(-abs(z), lower.tail = TRUE, log.p = TRUE)
>> >
>> > that could do really really well for large z compared to Rmpfr. Maybe I
>> am asking too much since
>> >
>> > z <-2
>> >> Rmpfr::format(2*pnorm(mpfr(-abs(z),100),lower.tail=TRUE,log.p=FALSE))
>> > [1] "1.660579603192917090365313727164e-86858901"
>> >
>> > already gives a rarely seen small p value. I gather I also need a
>> multiple precision exp() and their sum since exp(z^2/2) is also a Bayes
>> Factor so I get log(x_i )/sum_i log(x_i) instead. To this point, I am
>> obliged to clarify, see
>>
https://statgen.github.io/gwas-credible-sets/method/locuszoom-credible-sets.pdf
>> .
>> >
>> > I agree many feel geneticists go to far with small p values which I
>> would have difficulty to argue againston the other hand it is also
expected
>> to see these in a non-genetic context. For instance the Framingham study
>> was established in 1948 just got $34m for six years on phenotypewide
>> association which we would be interesting to see.
>> >
>> > Best wishes,
>> >
>> >
>> > Jing Hua
>> >
>> >
>> >
>> > From: peter dalgaard
>> > Sent: 21 June 2019 16:24
>> > To: jing hua zhao
>> > Cc: Rui Barradas; r-devel@r-project.org
>> > Subject: Re: [Rd] Calculation of e^{z^2/2} for a normal deviate z
>> >
>> > You may want to look into using the log option to qnorm
>> >
>> > e.g., in round figures:
>> >
>> >> log(1e-300)
>> > [1] -
Re: [Rd] Calculation of e^{z^2/2} for a normal deviate z
> jing hua zhao
> on Mon, 24 Jun 2019 08:51:43 + writes:
> Hi All,
> Thanks for all your comments which allows me to appreciate more of these
in Python and R.
> I just came across the matrixStats package,
> ## EXAMPLE #1
> lx <- c(1000.01, 1000.02)
> y0 <- log(sum(exp(lx)))
> print(y0) ## Inf
> y1 <- logSumExp(lx)
> print(y1) ## 1000.708
> and
>> ly <- lx*10
>> ly
> [1] 11000 12000
>> y1 <- logSumExp(ly)
>> print(y1)
> [1] 12000
>> logSumExp
> function (lx, idxs = NULL, na.rm = FALSE, ...)
> {
> has_na <- TRUE
> .Call(C_logSumExp, as.numeric(lx), idxs, as.logical(na.rm),
> has_na)
> }
>
>
> Maybe this is rather close?
Thank you Jing Hua,
indeed the issue of sums of (very large or very small)
exponentials is a special case, that can well be treated
specially
- as it is not so infrequent
- via "obvious" simplifications can be implemented even more accurately
We (authors of the R package 'copula') have had a need for
these as well, in likelihood computation for Archimedean
copulas, and
have efficient R level implementations, both for your case and
the -- even more delicate -- case of e.g., alternating signs of
exponential terms.
"Unfortunately", we had never exported the functions from the
package, so you'd need
copula:::lsum() # log sum {of exponentials in log scale}
or
copula:::lssum() # log *s*igned sum {of exponentials in log scale}
for the 2nd case.
The advantage is it's simple R code implemented quite
efficiently for the vector and matrices cases,
Source code from source file copula/R/special-func.R
(in svn at R-forge :
-->
https://r-forge.r-project.org/scm/viewvc.php/pkg/copula/R/special-func.R?view=markup&root=copula
)
{Yes, this GPL-2 licenced {with Copyright, see file, please
keep this one line}
## Copyright (C) 2012 Marius Hofert, Ivan Kojadinovic, Martin Maechler, and Jun
Y
}
--
##' Properly compute log(x_1 + .. + x_n) for a given (n x d)-matrix of n row
##' vectors log(x_1),..,log(x_n) (each of dimension d)
##' Here, x_i > 0 for all i
##' @title Properly compute the logarithm of a sum
##' @param lx (n,d)-matrix containing the row vectors log(x_1),..,log(x_n)
##'each of dimension d
##' @param l.off the offset to subtract and re-add; ideally in the order of
##'the maximum of each column
##' @return log(x_1 + .. + x_n) [i.e., OF DIMENSION d!!!] computed via
##' log(sum(x)) = log(sum(exp(log(x
##' = log(exp(log(x_max))*sum(exp(log(x)-log(x_max
##' = log(x_max) + log(sum(exp(log(x)-log(x_max)
##' = lx.max + log(sum(exp(lx-lx.max)))
##' => VECTOR OF DIMENSION d
##' @author Marius Hofert, Martin Maechler
lsum <- function(lx, l.off) {
rx <- length(d <- dim(lx))
if(mis.off <- missing(l.off)) l.off <- {
if(rx <= 1L)
max(lx)
else if(rx == 2L)
apply(lx, 2L, max)
}
if(rx <= 1L) { ## vector
if(is.finite(l.off))
l.off + log(sum(exp(lx - l.off)))
else if(mis.off || is.na(l.off) || l.off == max(lx))
l.off # NA || NaN or all lx == -Inf, or max(.) == Inf
else
stop("'l.off is infinite but not == max(.)")
} else if(rx == 2L) { ## matrix
if(any(x.off <- !is.finite(l.off))) {
if(mis.off || isTRUE(all.equal(l.off, apply(lx, 2L, max {
## we know l.off = colMax(.)
if(all(x.off)) return(l.off)
r <- l.off
iok <- which(!x.off)
l.of <- l.off[iok]
r[iok] <- l.of + log(colSums(exp(lx[,iok,drop=FALSE] -
rep(l.of, each=d[1]
r
} else ## explicitly specified l.off differing from colMax(.)
stop("'l.off' has non-finite values but differs from default
max(.)")
}
else
l.off + log(colSums(exp(lx - rep(l.off, each=d[1]
} else stop("not yet implemented for arrays of rank >= 3")
}
##' Properly compute log(x_1 + .. + x_n) for a given matrix of column vectors
##' log(|x_1|),.., log(|x_n|) and corresponding signs sign(x_1),.., sign(x_n)
##' Here, x_i is of arbitrary sign
##' @title compute logarithm of a sum with signed large coefficients
##' @param lxabs (d,n)-matrix containing the column vectors
log(|x_1|),..,log(|x_n|)
##'each of dimension d
##' @param signs corresponding matrix of signs sign(x_1), .., sign(x_n)
##' @param l.off the offset to subtract and re-add; ideally in the order of
max(.)
##' @param strict logical indicating if it should stop on some negative sums
##' @return log(x_1 + .. + x_n) [i.e., of dimension d] computed via
##' log(sum(x)) = log(sum(sign(x)*|x|)) =
log(sum(sign(x)*
Re: [Rd] Calculation of e^{z^2/2} for a normal deviate z
Hi Martin,
Thanks for another tips -- it will be hard for me to comprehend fully your
implementation of Rmpfr. I will look at these instead. I realised what I
thought was a difficult problem turned to be a popular one. I managed to get a
Python version here for information.
# http://bayesjumping.net/log-sum-exp-trick/
import numpy as np
def logSumExp(ns):
max = np.max(ns)
ds = ns - max
sumOfExp = np.exp(ds).sum()
return max + np.log(sumOfExp)
x = [11000, 12000]
logSumExp(x)
from scipy.misc import logsumexp
logsumexp(x)
import numpy as np;
prob1 = np.log(1e-50)
prob2 = np.log(2.5e-50)
prob12 = np.logaddexp(prob1, prob2)
prob12
np.exp(prob12)
As expected they are in good agreement with R.
Best regards,
Jing Hua
From: Martin Maechler
Sent: 24 June 2019 09:29
To: jing hua zhao
Cc: William Dunlap; Martin Maechler; r-devel@r-project.org
Subject: Re: [Rd] Calculation of e^{z^2/2} for a normal deviate z
> jing hua zhao
> on Mon, 24 Jun 2019 08:51:43 + writes:
> Hi All,
> Thanks for all your comments which allows me to appreciate more of these
in Python and R.
> I just came across the matrixStats package,
> ## EXAMPLE #1
> lx <- c(1000.01, 1000.02)
> y0 <- log(sum(exp(lx)))
> print(y0) ## Inf
> y1 <- logSumExp(lx)
> print(y1) ## 1000.708
> and
>> ly <- lx*10
>> ly
> [1] 11000 12000
>> y1 <- logSumExp(ly)
>> print(y1)
> [1] 12000
>> logSumExp
> function (lx, idxs = NULL, na.rm = FALSE, ...)
> {
> has_na <- TRUE
> .Call(C_logSumExp, as.numeric(lx), idxs, as.logical(na.rm),
> has_na)
> }
>
>
> Maybe this is rather close?
Thank you Jing Hua,
indeed the issue of sums of (very large or very small)
exponentials is a special case, that can well be treated
specially
- as it is not so infrequent
- via "obvious" simplifications can be implemented even more accurately
We (authors of the R package 'copula') have had a need for
these as well, in likelihood computation for Archimedean
copulas, and
have efficient R level implementations, both for your case and
the -- even more delicate -- case of e.g., alternating signs of
exponential terms.
"Unfortunately", we had never exported the functions from the
package, so you'd need
copula:::lsum() # log sum {of exponentials in log scale}
or
copula:::lssum() # log *s*igned sum {of exponentials in log scale}
for the 2nd case.
The advantage is it's simple R code implemented quite
efficiently for the vector and matrices cases,
Source code from source file copula/R/special-func.R
(in svn at R-forge :
-->
https://r-forge.r-project.org/scm/viewvc.php/pkg/copula/R/special-func.R?view=markup&root=copula
)
{Yes, this GPL-2 licenced {with Copyright, see file, please
keep this one line}
## Copyright (C) 2012 Marius Hofert, Ivan Kojadinovic, Martin Maechler, and Jun
Y
}
--
##' Properly compute log(x_1 + .. + x_n) for a given (n x d)-matrix of n row
##' vectors log(x_1),..,log(x_n) (each of dimension d)
##' Here, x_i > 0 for all i
##' @title Properly compute the logarithm of a sum
##' @param lx (n,d)-matrix containing the row vectors log(x_1),..,log(x_n)
##'each of dimension d
##' @param l.off the offset to subtract and re-add; ideally in the order of
##'the maximum of each column
##' @return log(x_1 + .. + x_n) [i.e., OF DIMENSION d!!!] computed via
##' log(sum(x)) = log(sum(exp(log(x
##' = log(exp(log(x_max))*sum(exp(log(x)-log(x_max
##' = log(x_max) + log(sum(exp(log(x)-log(x_max)
##' = lx.max + log(sum(exp(lx-lx.max)))
##' => VECTOR OF DIMENSION d
##' @author Marius Hofert, Martin Maechler
lsum <- function(lx, l.off) {
rx <- length(d <- dim(lx))
if(mis.off <- missing(l.off)) l.off <- {
if(rx <= 1L)
max(lx)
else if(rx == 2L)
apply(lx, 2L, max)
}
if(rx <= 1L) { ## vector
if(is.finite(l.off))
l.off + log(sum(exp(lx - l.off)))
else if(mis.off || is.na(l.off) || l.off == max(lx))
l.off # NA || NaN or all lx == -Inf, or max(.) == Inf
else
stop("'l.off is infinite but not == max(.)")
} else if(rx == 2L) { ## matrix
if(any(x.off <- !is.finite(l.off))) {
if(mis.off || isTRUE(all.equal(l.off, apply(lx, 2L, max {
## we know l.off = colMax(.)
if(all(x.off)) return(l.off)
r <- l.off
iok <- which(!x.off)
l.of <- l.off[iok]
r[iok] <- l.of + log(colSums(exp(lx[,iok,drop=FALSE] -
rep(l.of, each=d[1]
r
} else ## explicitly specified l.off differing from colMax(.)
Re: [Rd] methods package: A _R_CHECK_LENGTH_1_LOGIC2_=true error
> Henrik Bengtsson via R-core
> on Sun, 23 Jun 2019 11:29:58 -0700 writes:
> Thank you.
> To correct myself, I can indeed reproduce this with R --vanilla too.
> A reproducible example is:
> $ R --vanilla
> R version 3.6.0 Patched (2019-05-31 r76629) -- "Planting of a Tree"
> ...
>> Sys.setenv("_R_CHECK_LENGTH_1_LOGIC2_" = "true")
>> loadNamespace("oligo")
> Error in omittedSig && (signature[omittedSig] != "missing") :
> 'length(x) = 4 > 1' in coercion to 'logical(1)'
> Error: unable to load R code in package ‘oligo’
> /Henrik
Thank you Henrik, for the report, etc, but
hmm... after loading the oligo package, almost 40 (non
base+Recommended) packages have been loaded as well, which hence
need to have been installed before, too ..
which is not quite a "vanilla repr.ex." in my view
Worse, I cannot reproduce :
> Sys.setenv("_R_CHECK_LENGTH_1_LOGIC2_" = "true")
> Sys.getenv("_R_CHECK_LENGTH_1_LOGIC2_")
[1] "true"
> debugonce(conformMethod)
> loadNamespace("oligo")
Warning messages:
1: multiple methods tables found for ‘rowSums’
2: multiple methods tables found for ‘colSums’
3: multiple methods tables found for ‘rowMeans’
4: multiple methods tables found for ‘colMeans’
> sessionInfo()
R Under development (unstable) (2019-06-20 r76729)
(similarly with other versions of R >= 3.6.0).
So, even though R core has fixed this now in the sources, it
would be nice to have an "as simple as possible" repr.ex. for this.
Martin
> On Sun, Jun 23, 2019 at 1:54 AM peter dalgaard wrote:
>>
>> This looks obvious enough, so I just committed your fix to R-devel and
R-patched.
>>
>> I'm at the wrong machine for thorough testing, but at least it seems to
build OK. However, I sense some risk that this could uncover sleeping bugs
elsewhere, so watch out.
>>
>> -pd
>>
>> > On 22 Jun 2019, at 18:49 , Henrik Bengtsson
wrote:
>> >
>> > DISCLAIMER: I can not get this error with R --vanilla, so it only
>> > occurs when some other package is also loaded. I don't have time to
>> > find to narrow that down for a reproducible example, but I believe the
>> > following error in R 3.6.0:
>> >
>> >> Sys.setenv("_R_CHECK_LENGTH_1_LOGIC2_" = "true")
>> >> library(oligo)
>> > Error in omittedSig && (signature[omittedSig] != "missing") :
>> > 'length(x) = 4 > 1' in coercion to 'logical(1)'
>> > Error: unable to load R code in package 'oligo'
>> >
>> > is because of a '_R_CHECK_LENGTH_1_LOGIC2_=true' mistake in the
>> > 'methods' package. Here's the patch:
>> >
>> > $ svn diff src/library/methods/R/RMethodUtils.R &
>> > [1] 1062
>> > Index: src/library/methods/R/RMethodUtils.R
>> > ===
>> > --- src/library/methods/R/RMethodUtils.R (revision 76731)
>> > +++ src/library/methods/R/RMethodUtils.R (working copy)
>> > @@ -343,7 +343,7 @@
>> > call. = TRUE, domain = NA)
>> > }
>> > else if(!all(signature[omittedSig] == "missing")) {
>> > -omittedSig <- omittedSig && (signature[omittedSig] !=
"missing")
>> > +omittedSig <- omittedSig & (signature[omittedSig] !=
"missing")
>> > .message("Note: ", .renderSignature(f, sig0),
>> > gettextf("expanding the signature to include omitted
>> > arguments in definition: %s",
>> > paste(sigNames[omittedSig], "=
>> > \"missing\"",collapse = ", ")))
>> > [1]+ Donesvn diff
src/library/methods/R/RMethodUtils.R
>> >
>> > Maybe still in time for R 3.6.1?
>> >
>> > /Henrik
>> >
>> > __
>> > R-devel@r-project.org mailing list
>> > https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>> --
>> Peter Dalgaard, Professor,
>> Center for Statistics, Copenhagen Business School
>> Solbjerg Plads 3, 2000 Frederiksberg, Denmark
>> Phone: (+45)38153501
>> Office: A 4.23
>> Email: pd@cbs.dk Priv: pda...@gmail.com
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] methods package: A _R_CHECK_LENGTH_1_LOGIC2_=true error
**Maybe this bug needs to be understood further before applying the
patch because patch is most likely also wrong**
Because, from just looking at the expressions, I think neither the R
3.6.0 version:
omittedSig <- omittedSig && (signature[omittedSig] != "missing")
nor the patched version (I proposed):
omittedSig <- omittedSig & (signature[omittedSig] != "missing")
can be correct. For a starter, 'omittedSig' is a logical vector. We
see that 'omittedSig' is used to subset 'signature'. In other words,
the length of 'signature[omittedSig]' and hence the length of
'(signature[omittedSig] != "missing")' will equal sum(omittedSig),
i.e. the length will be in {0,1,...,length(omittedSig)}.
The R 3.6.0 version with '&&' is not correct because '&&' requires
length(omittedSig) == 1L and sum(omittedSig) == 1L, which is unlikely
to be the original intention.
The patched version with '&' is most likely not correct either because
'LHS & RHS' assume length(LHS) == length(RHS), unless it relies on the
auto-expansion of either vector. So, for the '&' version to be
correct, it basically requires that length(omittedSig) = length(LHS) =
length(RHS) = sum(omittedSig), which also sounds unlikely to be the
original intention.
Disclaimer: Please note that I have not at all studied the rest of the
function, so the above is just based on that single line plus
debugging that 'omittedSig' is a logical vector.
Martin, I don't have the time to dive into this further. Though I did
try to see if it happened when one of oligo's dependencies were
loaded, but that was not the case. It kicks in when oligo is loaded.
FYI, I can also replicate your non-replicatation on another R 3.6.0
version. I'll see if I can narrow down what's different, e.g.
comparing sessionInfo():s, etc. However, I want to reply with
findings above asap due to the R 3.6.1 release and you might roll back
the patch (since it might introduce other bugs as Peter mentioned).
/Henrik
On Mon, Jun 24, 2019 at 3:05 AM Martin Maechler
wrote:
>
> > Henrik Bengtsson via R-core
> > on Sun, 23 Jun 2019 11:29:58 -0700 writes:
>
> > Thank you.
> > To correct myself, I can indeed reproduce this with R --vanilla too.
> > A reproducible example is:
>
> > $ R --vanilla
> > R version 3.6.0 Patched (2019-05-31 r76629) -- "Planting of a Tree"
> > ...
> >> Sys.setenv("_R_CHECK_LENGTH_1_LOGIC2_" = "true")
> >> loadNamespace("oligo")
> > Error in omittedSig && (signature[omittedSig] != "missing") :
> > 'length(x) = 4 > 1' in coercion to 'logical(1)'
> > Error: unable to load R code in package ‘oligo’
>
> > /Henrik
>
> Thank you Henrik, for the report, etc, but
> hmm... after loading the oligo package, almost 40 (non
> base+Recommended) packages have been loaded as well, which hence
> need to have been installed before, too ..
> which is not quite a "vanilla repr.ex." in my view
>
> Worse, I cannot reproduce :
>
> > Sys.setenv("_R_CHECK_LENGTH_1_LOGIC2_" = "true")
> > Sys.getenv("_R_CHECK_LENGTH_1_LOGIC2_")
> [1] "true"
> > debugonce(conformMethod)
> > loadNamespace("oligo")
>
> Warning messages:
> 1: multiple methods tables found for ‘rowSums’
> 2: multiple methods tables found for ‘colSums’
> 3: multiple methods tables found for ‘rowMeans’
> 4: multiple methods tables found for ‘colMeans’
> > sessionInfo()
> R Under development (unstable) (2019-06-20 r76729)
>
> (similarly with other versions of R >= 3.6.0).
>
> So, even though R core has fixed this now in the sources, it
> would be nice to have an "as simple as possible" repr.ex. for this.
>
> Martin
>
>
>
> > On Sun, Jun 23, 2019 at 1:54 AM peter dalgaard wrote:
> >>
> >> This looks obvious enough, so I just committed your fix to R-devel and
> R-patched.
> >>
> >> I'm at the wrong machine for thorough testing, but at least it seems
> to build OK. However, I sense some risk that this could uncover sleeping bugs
> elsewhere, so watch out.
> >>
> >> -pd
> >>
> >> > On 22 Jun 2019, at 18:49 , Henrik Bengtsson
> wrote:
> >> >
> >> > DISCLAIMER: I can not get this error with R --vanilla, so it only
> >> > occurs when some other package is also loaded. I don't have time to
> >> > find to narrow that down for a reproducible example, but I believe
> the
> >> > following error in R 3.6.0:
> >> >
> >> >> Sys.setenv("_R_CHECK_LENGTH_1_LOGIC2_" = "true")
> >> >> library(oligo)
> >> > Error in omittedSig && (signature[omittedSig] != "missing") :
> >> > 'length(x) = 4 > 1' in coercion to 'logical(1)'
> >> > Error: unable to load R code in package 'oligo'
> >> >
> >> > is because of a '_R_CHECK_LENGTH_1_LOGIC2_=true' mistake in the
> >> > 'methods' package. Here's the patch:
> >> >
> >> > $ svn diff src/library/methods/R/RMethodUtils.R &
> >> > [1] 1062
> >> > Index: src/library/methods/R/RMethodUtils.R
> >> >
