Re: Low hanging fruit - Accelerated random distribution functions

[David Bateman]

Daniel J Sebald wrote:
> David Bateman wrote:
>> Paul Kienzle wrote:
>>
>>>> lognrnd,
>>>
>>> Current code can be improved a bit from:
>>>
>>> rnd(k) = exp(mu(k)) .* exp(sigma(k).*randn(1,length(k)));
>>>
>>> to:
>>>
>>> rnd(k) = exp(mu(k) + sigma(k).*randn(1,length(k)));
>>>
>>
>> This will only be a win in the case where mu or sigma are a matrix.
>
> It's still a win even for scalar mu and sigma, isn't it?  A pair of
> exp and multiply loses to an add and exp.  The extra exp is
> significant, probably the overhead for the octave interpretation masks
> that.  Anyway...
But its an exp of a scalar plus a scalar multiply, which is
insignificant.. Which is why I couldn't get an advantage

>
> Why the restriction of mu > 0?  The underlying exponentially
> distributed r.v. could have a negative mean with no consequence. 
> Naturally, sigma needs to be greater than 0 from a definition
> standpoint.  (And I wonder about the > 0 part because if variance is
> zero, that simply means the variable is a constant, exp(mu)).  And
> then even if the user specifies a sigma less than zero, it is
> effectively still a log normal random variable.  I.e., flip the sign of
>
> sigma(k) .* randn (1, length(k))
>
> and it is an r.v. with the same distribution.  (OK, may still want to
> take out the sigma < 0 case.)
>
> Also, I wonder if the special conditioning is necessary.  If mu = Inf,
> exp(Inf) = Inf.  If mu = -Inf, exp(-Inf) = 0.  How about avoiding the
> use of the index?  In most cases I'd think that the user will have no
> pathological cases.  I.e.,
>
>     rnd = exp(mu + sigma .* randn (1, length(sigma)));     k = find
> ((sigma < 0) | (sigma == Inf));
>     if (any (k))
>     rnd(k) = NaN
>     endif
>
> instead of
>
>     k = find ((mu > 0) & (mu < Inf) & (sigma > 0) & (sigma < Inf));
>     if (any (k))
>       rnd(k) = exp(mu(k) + sigma(k) .* randn (1, length(k)));     endif
>
> Does that help any for when k is the full array?

The restriction are reproduced from what was in the original code for
geometric_rnd. Looking at what the matlab code does it restricts sigma >
1 and has no restrictions on mu.. So yes I think you are right and we
should remove the restriction.

Yes I  agree that some of the special casing for NaN values might be
simplified, especially if the underlying operators return NaN for these
cases. I didn't both doing this at this point as I wanted the gains and
simplifications with the minimum work.. I incorporated your ideas in the
attached version of the patch.. However, similar thinks should be down
throughout the *rnd.m functions (while respecting the NaN values), which
I don't propose to do at this time.. Want to propose a patch?

Regards
David



>
> Dan
>


-- 
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/lognrnd.m.orig46 2007-02-23 
11:46:17.949059912 +0100
--- ./scripts/statistics/distributions/lognrnd.m        2007-02-23 
14:10:35.399814087 +0100
***************
*** 78,100 ****
    endif
  
    if (isscalar (mu) && isscalar (sigma))
!     if  (!(mu > 0) | !(mu < Inf) | !(sigma > 0) | !(sigma < Inf))
        rnd = NaN * ones (sz);
-     elseif find ((mu > 0) & (mu < Inf) & (sigma > 0) & (sigma < Inf));
-       rnd = exp (mu) * exp (sigma .* randn (sz));
      else
!       rnd = zeros (sz);
      endif
    else
!     rnd = zeros (sz);
!     k = find (!(mu > 0) | !(mu < Inf) | !(sigma > 0) | !(sigma < Inf));
      if (any (k))
!       rnd(k) = NaN * ones (1, length (k));
!     endif
! 
!     k = find ((mu > 0) & (mu < Inf) & (sigma > 0) & (sigma < Inf));
!     if (any (k))
!       rnd(k) = exp (mu(k)) .* exp (sigma(k) .* randn (1, length (k)));
      endif
    endif
  
--- 78,93 ----
    endif
  
    if (isscalar (mu) && isscalar (sigma))
!     if  (!(sigma > 0) || !(sigma < Inf))
        rnd = NaN * ones (sz);
      else
!       rnd = exp(mu + sigma .* randn (sz)); 
      endif
    else
!     rnd = exp (mu + sigma .* randn (sz));
!     k = find ((sigma < 0) | (sigma == Inf));
      if (any (k))
!       rnd(k) = NaN
      endif
    endif