Yes, that looks like a bug and an easily fixable one too.
However, I spy another issue: Why do we check the !R_FINITE(x) && mu ==
x before checking for sd < 0 ? The difference is whether we
return ML_NAN;
or
ML_ERR_return_NAN;
but surely negative sd should always be an error?
I'd be inclined to do
if (sigma < 0) ML_ERR_return_NAN;
if(!R_FINITE(sigma)) return R_D__0;
if(!R_FINITE(x) && mu == x) return ML_NAN;/* x-mu is NaN */
if (sigma == 0)
return (x == mu) ? ML_POSINF : R_D__0;
x = (x - mu) / sigma;
(Ping Martin...)
-pd
> On 7 Dec 2019, at 23:40 , Wang Jiefei <szwjf08 at gmail.com> wrote:
>
> Good question, I cannot speak for R's developers but I would like to
> provide some information on the problem. Here are the first few lines of
> the dnorm function located at src\nmath\dnorm.c:
>
> ```
> double dnorm4(double x, double mu, double sigma, int give_log)
> {
> #ifdef IEEE_754
> if (ISNAN(x) || ISNAN(mu) || ISNAN(sigma))
> return x + mu + sigma;
> #endif
> if(!R_FINITE(sigma)) return R_D__0;
> if(!R_FINITE(x) && mu == x) return ML_NAN;/* x-mu is NaN */
> if (sigma <= 0) {
> if (sigma < 0) ML_ERR_return_NAN;
> /* sigma == 0 */
> return (x == mu) ? ML_POSINF : R_D__0;
> }
> ....
> }
> ```
>
> You can clearly see where the problem is. I think either the document or
> the code needs a modification.
>
> Best,
> Jiefei
>
> On Sat, Dec 7, 2019 at 5:05 PM Weigand, Stephen D. via R-devel <
> r-devel at r-project.org> wrote:
>
>> Hi,
>> Apropos of a recent Inf question, I've previously wondered if dnorm
"does
>> the right thing" with
>>
>> dnorm(0, 0, -Inf)
>>
>> which gives zero. Should that be zero or NaN (or NA)?
>>
>> The help says "'sd < 0' is an error and returns
'NaN'" and since -Inf < 0
>> is TRUE, then... is this a bug?
>>
>> Thank you,
>> Stephen
>> Rochester, MN USA
>>
>> ______________________________________________
>> R-devel at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-devel at 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.mes at cbs.dk Priv: PDalgd at gmail.com
>>>>> peter dalgaard >>>>> on Sun, 8 Dec 2019 12:11:50 +0100 writes:> Yes, that looks like a bug and an easily fixable one too. agreed. > However, I spy another issue: Why do we check the > !R_FINITE(x) && mu == x before checking for sd < 0 ? The > difference is whether we > return ML_NAN; or ML_ERR_return_NAN; > but surely negative sd should always be an error? > I'd be inclined to do > if (sigma < 0) ML_ERR_return_NAN; > if(!R_FINITE(sigma)) return R_D__0; > if(!R_FINITE(x) && mu == x) return ML_NAN;/* x-mu is NaN */ > if (sigma == 0) > return (x == mu) ? ML_POSINF : R_D__0; > x = (x - mu) / sigma; > (Ping Martin...) I think you are spot on, Peter. All of this code has a longish history, with incremental border case improvements. Let's hope (somewhat unrealistically) this is the last one for dnorm(). NB: dlnorm() and some of the gamma/chisq/.. may need a similar adjustment Lastly, regression tests for this (either in tests/d-p-q-r-tests.{R,Rout.save} or easier in reg-tests-1d.R) should be added too. > -pd >> On 7 Dec 2019, at 23:40 , Wang Jiefei <szwjf08 at gmail.com> wrote: >> >> Good question, I cannot speak for R's developers but I would like to >> provide some information on the problem. Here are the first few lines of >> the dnorm function located at src\nmath\dnorm.c: >> >> ``` >> double dnorm4(double x, double mu, double sigma, int give_log) >> { >> #ifdef IEEE_754 >> if (ISNAN(x) || ISNAN(mu) || ISNAN(sigma)) >> return x + mu + sigma; >> #endif >> if(!R_FINITE(sigma)) return R_D__0; >> if(!R_FINITE(x) && mu == x) return ML_NAN;/* x-mu is NaN */ >> if (sigma <= 0) { >> if (sigma < 0) ML_ERR_return_NAN; >> /* sigma == 0 */ >> return (x == mu) ? ML_POSINF : R_D__0; >> } >> .... >> } >> ``` >> >> You can clearly see where the problem is. I think either the document or >> the code needs a modification. >> >> Best, >> Jiefei >> >> On Sat, Dec 7, 2019 at 5:05 PM Weigand, Stephen D. via R-devel < >> r-devel at r-project.org> wrote: >> >>> Hi, >>> Apropos of a recent Inf question, I've previously wondered if dnorm "does >>> the right thing" with >>> >>> dnorm(0, 0, -Inf) >>> >>> which gives zero. Should that be zero or NaN (or NA)? >>> >>> The help says "'sd < 0' is an error and returns 'NaN'" and since -Inf < 0 >>> is TRUE, then... is this a bug? >>> >>> Thank you, >>> Stephen >>> Rochester, MN USA >>> >>> ______________________________________________ >>> R-devel at r-project.org mailing list >>> https://stat.ethz.ch/mailman/listinfo/r-devel >>> >> >> [[alternative HTML version deleted]] >> >> ______________________________________________ >> R-devel at 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.mes at cbs.dk Priv: PDalgd at gmail.com
I have committed a fix for r-devel (dnorm only). -pd> On 9 Dec 2019, at 08:49 , Martin Maechler <maechler at stat.math.ethz.ch> wrote: > >>>>>> peter dalgaard >>>>>> on Sun, 8 Dec 2019 12:11:50 +0100 writes: > >> Yes, that looks like a bug and an easily fixable one too. > > agreed. > >> However, I spy another issue: Why do we check the >> !R_FINITE(x) && mu == x before checking for sd < 0 ? The >> difference is whether we > >> return ML_NAN; or ML_ERR_return_NAN; > >> but surely negative sd should always be an error? > >> I'd be inclined to do > >> if (sigma < 0) ML_ERR_return_NAN; >> if(!R_FINITE(sigma)) return R_D__0; >> if(!R_FINITE(x) && mu == x) return ML_NAN;/* x-mu is NaN */ >> if (sigma == 0) >> return (x == mu) ? ML_POSINF : R_D__0; >> x = (x - mu) / sigma; > > >> (Ping Martin...) > > I think you are spot on, Peter. > All of this code has a longish history, with incremental border > case improvements. > Let's hope (somewhat unrealistically) this is the last one for > dnorm(). > > NB: dlnorm() and some of the gamma/chisq/.. may need a > similar adjustment > > Lastly, regression tests for this > (either in tests/d-p-q-r-tests.{R,Rout.save} > or easier in reg-tests-1d.R) should be added too. > >> -pd > >>> On 7 Dec 2019, at 23:40 , Wang Jiefei <szwjf08 at gmail.com> wrote: >>> >>> Good question, I cannot speak for R's developers but I would like to >>> provide some information on the problem. Here are the first few lines of >>> the dnorm function located at src\nmath\dnorm.c: >>> >>> ``` >>> double dnorm4(double x, double mu, double sigma, int give_log) >>> { >>> #ifdef IEEE_754 >>> if (ISNAN(x) || ISNAN(mu) || ISNAN(sigma)) >>> return x + mu + sigma; >>> #endif >>> if(!R_FINITE(sigma)) return R_D__0; >>> if(!R_FINITE(x) && mu == x) return ML_NAN;/* x-mu is NaN */ >>> if (sigma <= 0) { >>> if (sigma < 0) ML_ERR_return_NAN; >>> /* sigma == 0 */ >>> return (x == mu) ? ML_POSINF : R_D__0; >>> } >>> .... >>> } >>> ``` >>> >>> You can clearly see where the problem is. I think either the document or >>> the code needs a modification. >>> >>> Best, >>> Jiefei >>> >>> On Sat, Dec 7, 2019 at 5:05 PM Weigand, Stephen D. via R-devel < >>> r-devel at r-project.org> wrote: >>> >>>> Hi, >>>> Apropos of a recent Inf question, I've previously wondered if dnorm "does >>>> the right thing" with >>>> >>>> dnorm(0, 0, -Inf) >>>> >>>> which gives zero. Should that be zero or NaN (or NA)? >>>> >>>> The help says "'sd < 0' is an error and returns 'NaN'" and since -Inf < 0 >>>> is TRUE, then... is this a bug? >>>> >>>> Thank you, >>>> Stephen >>>> Rochester, MN USA >>>> >>>> ______________________________________________ >>>> R-devel at r-project.org mailing list >>>> https://stat.ethz.ch/mailman/listinfo/r-devel >>>> >>> >>> [[alternative HTML version deleted]] >>> >>> ______________________________________________ >>> R-devel at 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.mes at cbs.dk Priv: PDalgd at gmail.com-- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Office: A 4.23 Email: pd.mes at cbs.dk Priv: PDalgd at gmail.com