Re: Overhaul of statistical distribution functions

[Ben Abbott]

On Sep 20, 2011, at 3:13 PM, Rik wrote:

> 9/20/11
> 
> All,
> 
> As another concrete achievement for the future 3.6 release, I've rewritten
> 89 statistical distribution functions in scripts/statistics/distributions. 
> Quoting directly from the NEWS file
> 
> ** Overhaul of statistical distribution functions
> 
>    Functions now return "single" outputs for inputs of class "single".
> 
>    75% reduction in memory usage through use of logical indexing.
> 
>    Random sample functions now use the same syntax as rand() and accept
>    a comma separated list of dimensions or a dimension vector.
> 
>    Functions have been made Matlab-compatible with regard to special
>    cases (probability on boundaries, probabilities for values outside
>    distribution, etc.).  This may cause subtle changes to existing
>    scripts.
> 
>    negative binomial function has been extended to real, non-integer inputs.
>    discrete_inv() now returns v(1) for 0 instead of NaN.
>    nbincdf() recoded to use closed form solution with betainc().
> 
> All of the rewritten functions have complete tests as well.  Thanks in
> particular to Michael Godfrey who has run dozens of odd corner cases in
> Matlab for me.  I'm checking in the changes now in the hope that some of
> the inevitable bugs this large a change provokes will be caught.
> 
> Cheers,
> Rik

Rik,

I'm happy to see you complete this effort.

With this change, I see 5 failures. They all appear to be due to tolerances. 
All tests pass if I use a tolerance of eps(). My diff is also below.

Is adding tolerance appropriate? If so, I'm happy to push a changeset.

Ben

>>>>> processing 
>>>>> /Users/bpabbott/Development/mercurial/default/scripts/statistics/distributions/binopdf.m
  ***** assert(binopdf (x, 2*ones(1,5), 0.5*ones(1,5)), y);
!!!!! test failed
assert (binopdf (x, 2 * ones (1, 5), 0.5 * ones (1, 5)),y) expected
   0.00000   0.25000   0.50000   0.25000   0.00000
but got
   0.00000   0.25000   0.50000   0.25000   0.00000
values do not match
shared variables   scalar structure containing the fields:
    x =
      -1   0   1   2   3
    y =
       0.00000   0.25000   0.50000   0.25000   0.00000
  ***** assert(binopdf (x, 2, 0.5*ones(1,5)), y);
!!!!! test failed
assert (binopdf (x, 2, 0.5 * ones (1, 5)),y) expected
   0.00000   0.25000   0.50000   0.25000   0.00000
but got
   0.00000   0.25000   0.50000   0.25000   0.00000
values do not match
shared variables   scalar structure containing the fields:
    x =
      -1   0   1   2   3
    y =
       0.00000   0.25000   0.50000   0.25000   0.00000
  ***** assert(binopdf (x, 2*ones(1,5), 0.5), y);
!!!!! test failed
assert (binopdf (x, 2 * ones (1, 5), 0.5),y) expected
   0.00000   0.25000   0.50000   0.25000   0.00000
but got
   0.00000   0.25000   0.50000   0.25000   0.00000
values do not match
shared variables   scalar structure containing the fields:
    x =
      -1   0   1   2   3
    y =
       0.00000   0.25000   0.50000   0.25000   0.00000
  ***** assert(binopdf ([x, NaN], 2, 0.5), [y, NaN]);
!!!!! test failed
assert (binopdf ([x, NaN], 2, 0.5),[y, NaN]) expected
   0.00000   0.25000   0.50000   0.25000   0.00000       NaN
but got
   0.00000   0.25000   0.50000   0.25000   0.00000       NaN
values do not match
shared variables   scalar structure containing the fields:
    x =
      -1   0   1   2   3
    y =
       0.00000   0.25000   0.50000   0.25000   0.00000


>>>>> processing 
>>>>> /Users/bpabbott/Development/mercurial/default/scripts/statistics/distributions/tpdf.m
  ***** assert(tpdf ([x, NaN], 1), [y, NaN]);
!!!!! test failed
assert (tpdf ([x, NaN], 1),[y, NaN]) expected
   0.00000   0.31831   0.25465   0.15915   0.00000       NaN
but got
   0.00000   0.31831   0.25465   0.15915   0.00000       NaN
values do not match
shared variables   scalar structure containing the fields:
    x =
          -Inf   0.00000   0.50000   1.00000       Inf
    y =
       0.00000   0.31831   0.25465   0.15915   0.00000

diff --git a/scripts/statistics/distributions/binopdf.m 
b/scripts/statistics/distributions/binopdf.m
--- a/scripts/statistics/distributions/binopdf.m
+++ b/scripts/statistics/distributions/binopdf.m
@@ -66,15 +66,16 @@
 endfunction
 
 
-%!shared x,y
+%!shared x, y, tol
 %! x = [-1 0 1 2 3];
 %! y = [0 1/4 1/2 1/4 0];
-%!assert(binopdf (x, 2*ones(1,5), 0.5*ones(1,5)), y);
-%!assert(binopdf (x, 2, 0.5*ones(1,5)), y);
-%!assert(binopdf (x, 2*ones(1,5), 0.5), y);
+%! tol = eps;
+%!assert(binopdf (x, 2*ones(1,5), 0.5*ones(1,5)), y, tol);
+%!assert(binopdf (x, 2, 0.5*ones(1,5)), y, tol);
+%!assert(binopdf (x, 2*ones(1,5), 0.5), y, tol);
 %!assert(binopdf (x, 2*[0 -1 NaN 1.1 1], 0.5), [0 NaN NaN NaN 0]);
 %!assert(binopdf (x, 2, 0.5*[0 -1 NaN 3 1]), [0 NaN NaN NaN 0]);
-%!assert(binopdf ([x, NaN], 2, 0.5), [y, NaN]);
+%!assert(binopdf ([x, NaN], 2, 0.5), [y, NaN], tol);
 
 %% Test class of input preserved
 %!assert(binopdf (single([x, NaN]), 2, 0.5), single([y, NaN]));
diff --git a/scripts/statistics/distributions/tpdf.m 
b/scripts/statistics/distributions/tpdf.m
--- a/scripts/statistics/distributions/tpdf.m
+++ b/scripts/statistics/distributions/tpdf.m
@@ -78,7 +78,7 @@
 %!assert(tpdf (x, [0 NaN 1 1 1]), [NaN NaN y(3:5)], eps);
 
 %% Test class of input preserved
-%!assert(tpdf ([x, NaN], 1), [y, NaN]);
+%!assert(tpdf ([x, NaN], 1), [y, NaN], eps);
 %!assert(tpdf (single([x, NaN]), 1), single([y, NaN]), eps("single"));
 %!assert(tpdf ([x, NaN], single(1)), single([y, NaN]), eps("single"));