What should dnorm(0, 0, -Inf) return?

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

What should dnorm(0, 0, -Inf) return?

R devel mailing list
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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: What should dnorm(0, 0, -Inf) return?

Wang Jiefei
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 <
[hidden email]> 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
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: What should dnorm(0, 0, -Inf) return?

Peter Dalgaard-2
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 <[hidden email]> 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 <
> [hidden email]> 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
>>
>> ______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] 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: [hidden email]  Priv: [hidden email]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: What should dnorm(0, 0, -Inf) return?

Martin Maechler
>>>>> 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 <[hidden email]> 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 <
    >> [hidden email]> 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
    >>>
    >>> ______________________________________________
    >>> [hidden email] mailing list
    >>> https://stat.ethz.ch/mailman/listinfo/r-devel
    >>>
    >>
    >> [[alternative HTML version deleted]]
    >>
    >> ______________________________________________
    >> [hidden email] 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: [hidden email]  Priv: [hidden email]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: What should dnorm(0, 0, -Inf) return?

Peter Dalgaard-2
I have committed a fix for r-devel (dnorm only).

-pd

> On 9 Dec 2019, at 08:49 , Martin Maechler <[hidden email]> 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 <[hidden email]> 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 <
>>> [hidden email]> 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
>>>>
>>>> ______________________________________________
>>>> [hidden email] mailing list
>>>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>>>
>>>
>>> [[alternative HTML version deleted]]
>>>
>>> ______________________________________________
>>> [hidden email] 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: [hidden email]  Priv: [hidden email]

--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: [hidden email]  Priv: [hidden email]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel