[Octave-bug-tracker] [bug #40470] Wrong output from normcdf (with exp?) for negatives

[anonymous]

URL:
  <http://savannah.gnu.org/bugs/?40470>

                 Summary: Wrong output from normcdf (with exp?) for  negatives

                 Project: GNU Octave
            Submitted by: None
            Submitted on: Mon 04 Nov 2013 03:04:56 PM UTC
                Category: Octave Function
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: Inaccurate Result
                  Status: None
             Assigned to: None
         Originator Name: Alan Genz
        Originator Email: address@hidden
             Open/Closed: Open
         Discussion Lock: Any
                 Release: 3.2.4
        Operating System: GNU/Linux

    _______________________________________________________

Details:

octave:21> for t=-[7:10],disp([t exp(-t^2/2) normcdf(t)]),end
  -7.0000e+00   2.2897e-11   1.2798e-12
  -8.0000e+00   1.2664e-14   6.1062e-16
  -9.00000   0.00000   0.00000
  -10.00000    0.00000    0.00000
octave:22> for t=-[7:10],disp([t exp(-t^2/2)]),end
  -7.0000e+00   2.2897e-11
  -8.0000e+00   1.2664e-14
  -9.0000e+00   2.5768e-18
  -1.0000e+01   1.9287e-22
octave:23> for t=-[7:10],disp([t normcdf(t)]),end
  -7.0000e+00   1.2798e-12
  -8.0000e+00   6.1062e-16
  -9   0
  -10    0

The above results are inconsistent; normcdf should not 
underflow to 0 until t is near -37. 
But why is exp also affected? 
This bug has been reported previously (see #29543).
The simple fix is to use normcdf(t) = erfc-t/sqrt(2))/2:
octave:27> for t=-[7:10],disp([t exp(-t^2/2) ...
                   erfc(-t/sqrt(2))/2]),end
  -7.0000e+00   2.2897e-11   1.2798e-12
  -8.0000e+00   1.2664e-14   6.2210e-16
  -9.0000e+00   2.5768e-18   1.1286e-19
  -1.0000e+01   1.9287e-22   7.6199e-24

Here are some Matlab results for comparison
>> for t=-[7:10],disp([t exp(-t^2/2) ...
     normcdf(t)]),end                     
           -7   2.2897e-11   1.2798e-12
           -8   1.2664e-14    6.221e-16
           -9   2.5768e-18   1.1286e-19
          -10   1.9287e-22   7.6199e-24






    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?40470>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/