Re: consistent crash using gammrnd

[David Bateman]

Hermann wrote:
> Am Samstag, 14. Juli 2007 schrieb Paul Kienzle:
>
>   
>>> gamrnd(rand(100,1),rand(100,1))
>>>       
>> The attached patch fixes the gamrnd bug ... it was doing
>> vector/vector rather than vector./vector.
>>     
>
> Hi,
> I’m not sure, but maybe this should be .* instead of ./ ? I noted that 
> in 'help randg' it sais
>
> | This can be used to generate many distributions:
> |
> | `gamma (a, b)' for `a > -1', `b > 0'
> |            r = b * randg (a)
>
>   
Well, yes and no.. The current definition of gampdf, gamcdf and gamrnd
use a definition of 1/mu for the parameter B. So the current behavior of
gamrnd is consistent with gamcdf, gampdf and gaminv. Unfortunately this
seems to be the reverse of what is done in Matlab. For compatibility its
probably better to mimic the matlab behavior. Consider the attached
patch that changes the behavior to be consistent with matlab.. Note that
exppdf, etc are also concerned are they are really just special cases of
the gamma distributions.

Do the other statistics functions have the same issues? I don't have the
statistics toolbox for matlab and so can't easily check...

D.

-- 
David Bateman                                address@hidden
Motorola Labs - Paris                        +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 6 72 01 06 33 (Mob) 
91193 Gif-Sur-Yvette FRANCE                  +33 1 69 35 77 01 (Fax) 

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary

*** ./scripts/statistics/distributions/gamrnd.m.orig43  2007-07-16 
11:17:08.208590511 +0200
--- ./scripts/statistics/distributions/gamrnd.m 2007-07-16 11:27:14.029964245 
+0200
***************
*** 82,88 ****
      if (find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf)))
        rnd = NaN * ones (sz);
      else
!       rnd = randg(a,sz)/b;
      endif
    else 
      k = find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf));
--- 82,88 ----
      if (find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf)))
        rnd = NaN * ones (sz);
      else
!       rnd = b .* randg(a, sz);
      endif
    else 
      k = find (!(a > 0) | !(a < Inf) | !(b > 0) | !(b < Inf));
***************
*** 91,97 ****
      endif
      k = find ((a > 0) & (a < Inf) & (b > 0) & (b < Inf));
      if (any (k))
!       rnd(k) = randg(a(k),size(k))/b(k);
      endif
    endif
  
--- 91,97 ----
      endif
      k = find ((a > 0) & (a < Inf) & (b > 0) & (b < Inf));
      if (any (k))
!       rnd(k) = b(k) .* randg(a(k), size(k));
      endif
    endif
  
*** ./scripts/statistics/distributions/gampdf.m.orig43  2007-07-16 
11:17:14.374260939 +0200
--- ./scripts/statistics/distributions/gampdf.m 2007-07-16 11:34:09.546037770 
+0200
***************
*** 52,73 ****
    k = find ((x > 0) & (a > 0) & (a <= 1) & (b > 0));
    if (any (k))
      if (isscalar(a) && isscalar(b))
!       pdf(k) = ((b .^ a) .* (x(k) .^ (a - 1))
!               .* exp(-b .* x(k)) ./ gamma (a));
      else
!       pdf(k) = ((b(k) .^ a(k)) .* (x(k) .^ (a(k) - 1))
!               .* exp(-b(k) .* x(k)) ./ gamma (a(k)));
      endif
    endif
  
    k = find ((x > 0) & (a > 1) & (b > 0));
    if (any (k))
      if (isscalar(a) && isscalar(b))
!       pdf(k) = exp (a .* log (b) + (a-1) .* log (x(k))
!                   - b .* x(k) - gammaln (a));
      else
!       pdf(k) = exp (a(k) .* log (b(k)) + (a(k)-1) .* log (x(k))
!                   - b(k) .* x(k) - gammaln (a(k)));
      endif
    endif
  
--- 52,73 ----
    k = find ((x > 0) & (a > 0) & (a <= 1) & (b > 0));
    if (any (k))
      if (isscalar(a) && isscalar(b))
!       pdf(k) = (x(k) .^ (a - 1)) ...
!               .* exp(- x(k) ./ b) ./ gamma (a) ./ (b .^ a);
      else
!       pdf(k) = (x(k) .^ (a(k) - 1)) ...
!               .* exp(- x(k) ./ b(k)) ./ gamma (a(k)) ./ (b(k) .^ a(k));
      endif
    endif
  
    k = find ((x > 0) & (a > 1) & (b > 0));
    if (any (k))
      if (isscalar(a) && isscalar(b))
!       pdf(k) = exp (- a .* log (b) + (a-1) .* log (x(k))
!                   - x(k) ./ b - gammaln (a));
      else
!       pdf(k) = exp (- a(k) .* log (b(k)) + (a(k)-1) .* log (x(k))
!                   - x(k) ./ b(k) - gammaln (a(k)));
      endif
    endif
  
*** ./scripts/statistics/distributions/gamcdf.m.orig43  2007-07-16 
11:17:17.798077943 +0200
--- ./scripts/statistics/distributions/gamcdf.m 2007-07-16 11:35:08.970067259 
+0200
***************
*** 52,60 ****
    k = find ((x > 0) & (a > 0) & (b > 0));
    if (any (k))
      if (isscalar (a) && isscalar(b))
!       cdf (k) = gammainc (b * x(k), a);
      else
!       cdf (k) = gammainc (b(k) .* x(k), a(k));
      endif
    endif
  
--- 52,60 ----
    k = find ((x > 0) & (a > 0) & (b > 0));
    if (any (k))
      if (isscalar (a) && isscalar(b))
!       cdf (k) = gammainc (x(k) ./ b, a);
      else
!       cdf (k) = gammainc (x(k) ./ b(k), a(k));
      endif
    endif
  
*** ./scripts/statistics/distributions/gaminv.m.orig43  2007-07-16 
11:17:37.781010129 +0200
--- ./scripts/statistics/distributions/gaminv.m 2007-07-16 11:36:39.753521593 
+0200
***************
*** 59,67 ****
      if (!isscalar(a) || !isscalar(b))
        a = a (k);
        b = b (k);
!       y = a ./ b;
      else
!       y = a / b * ones (size (k));
      endif
      x = x (k);
      l = find (x < eps);
--- 59,67 ----
      if (!isscalar(a) || !isscalar(b))
        a = a (k);
        b = b (k);
!       y = a .* b;
      else
!       y = a * b * ones (size (k));
      endif
      x = x (k);
      l = find (x < eps);
*** ./scripts/statistics/distributions/exprnd.m.orig43  2007-07-16 
11:19:23.840349088 +0200
--- ./scripts/statistics/distributions/exprnd.m 2007-07-16 11:19:09.946090120 
+0200
***************
*** 69,75 ****
  
    if (isscalar (l))
      if ((l > 0) && (l < Inf))
!       rnd = rande(sz) / l;
      else
        rnd = NaN * ones (sz);
      endif
--- 69,75 ----
  
    if (isscalar (l))
      if ((l > 0) && (l < Inf))
!       rnd = rande(sz) * l;
      else
        rnd = NaN * ones (sz);
      endif
***************
*** 81,87 ****
      endif
      k = find ((l > 0) & (l < Inf));
      if (any (k))
!       rnd(k) = rande(size(k)) / l(k);
      endif
    endif
  
--- 81,87 ----
      endif
      k = find ((l > 0) & (l < Inf));
      if (any (k))
!       rnd(k) = rande(size(k)) .* l(k);
      endif
    endif
  
*** ./scripts/statistics/distributions/exppdf.m.orig43  2007-07-16 
11:19:33.938810606 +0200
--- ./scripts/statistics/distributions/exppdf.m 2007-07-16 11:20:57.227414605 
+0200
***************
*** 54,64 ****
    k = find ((x > 0) & (x < Inf) & (l > 0));
    if (any (k))
      if isscalar (l)
!       pdf(k) = l .* exp (- l .* x(k));
      elseif isscalar (x)
!       pdf(k) = l(k) .* exp (- l(k) .* x);
      else
!       pdf(k) = l(k) .* exp (- l(k) .* x(k));
      endif
    endif
  
--- 54,64 ----
    k = find ((x > 0) & (x < Inf) & (l > 0));
    if (any (k))
      if isscalar (l)
!       pdf(k) = exp (- x(k) ./ l) ./ l;
      elseif isscalar (x)
!       pdf(k) = exp (- x ./ l(k)) ./ l(k);
      else
!       pdf(k) = exp (- x(k) ./ l(k)) ./ l(k);
      endif
    endif
  
*** ./scripts/statistics/distributions/expcdf.m.orig43  2007-07-16 
11:19:37.124640746 +0200
--- ./scripts/statistics/distributions/expcdf.m 2007-07-16 11:23:32.623406482 
+0200
***************
*** 63,73 ****
    k = find ((x > 0) & (x < Inf) & (l > 0));
    if (any (k))
      if isscalar (l)
!       cdf (k) = 1 - exp (- l .* x(k));
      elseif isscalar (x)
!       cdf (k) = 1 - exp (- l(k) .* x);
      else
!       cdf (k) = 1 - exp (- l(k) .* x(k));
      endif
    endif
  
--- 63,73 ----
    k = find ((x > 0) & (x < Inf) & (l > 0));
    if (any (k))
      if isscalar (l)
!       cdf (k) = 1 - exp (- x(k) ./ l);
      elseif isscalar (x)
!       cdf (k) = 1 - exp (- x ./ l(k));
      else
!       cdf (k) = 1 - exp (- x(k) ./ l(k));
      endif
    endif
  
*** ./scripts/statistics/distributions/expinv.m.orig43  2007-07-16 
11:19:42.978328665 +0200
--- ./scripts/statistics/distributions/expinv.m 2007-07-16 11:26:28.561316833 
+0200
***************
*** 61,71 ****
    k = find ((x > 0) & (x < 1) & (l > 0));
    if (any (k))
      if isscalar (l)
!       inv(k) = - log (1 - x(k)) ./ l;
      elseif isscalar (x)
!       inv(k) = - log (1 - x) ./ l(k);
      else
!       inv(k) = - log (1 - x(k)) ./ l(k);
      endif
    endif
  
--- 61,71 ----
    k = find ((x > 0) & (x < 1) & (l > 0));
    if (any (k))
      if isscalar (l)
!       inv(k) = - l .* log (1 - x(k));
      elseif isscalar (x)
!       inv(k) = - l(k) .* log (1 - x);
      else
!       inv(k) = - l(k) .* log (1 - x(k));
      endif
    endif
  
*** ./NEWS.orig43       2007-07-16 11:25:24.807613176 +0200
--- ./NEWS      2007-07-16 11:41:03.849249381 +0200
***************
*** 147,150 ****
--- 147,155 ----
      normrnd have been modified to compute distributions using the
      standard deviation instead of the variance.
  
+  ** For compatibility with Matlab, gamcdf, gaminv, gampdf, gamrnd,
+     expcdf, expinv, exppdf and exprnd have been modified to compute 
+     the distributions using the standard scale factor rather than
+     one over the scale factor.
+ 
  See NEWS.2 for old news.

2007-07-16  David Bateman  <address@hidden>

        * statistics/distributions/gamcdf.m, statistics/distributions/gaminv.m,
        statistics/distributions/gampdf.m, statistics/distributions/gamrnd.m,
        statistics/distributions/expcdf.m, statistics/distributions/expinv.m,
        statistics/distributions/exppdf.m, statistics/distributions/exprnd.m:
        Use standard scale factor rather than one on the scale factor for
        compatibility.