From pd@|gd @end|ng |rom gm@||@com  Mon Dec  9 12:59:07 2019
From: pd@|gd @end|ng |rom gm@||@com (peter dalgaard)
Date: Mon, 9 Dec 2019 12:59:07 +0100
Subject: [Rd] What should dnorm(0, 0, -Inf) return?
In-Reply-To: <24045.64646.983266.213654@stat.math.ethz.ch>
References: <28fddd$csrn33@ironport10.mayo.edu>
 <CAGiFhPP4AVZnZUfFo=N5CPJZCUSgL6sWUROaT8fXAFRvKqpFVA@mail.gmail.com>
 <2E7FB3D3-4FB2-49CF-983D-77735A79D7AD@gmail.com>
 <24045.64646.983266.213654@stat.math.ethz.ch>
Message-ID: <75FD9E0E-2BC1-493E-A4EF-2FE221D2DDAF@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



[Rd] What should dnorm(0, 0, -Inf) return?
peter dalgaard