This is not a big deal to me, but what exactly would break by allowing λ=0? Returning inf for skewness and kurtosis and an error for entropy doesn't seem problematic to me. I don't do estimation of KL divergences, but it appears a bit strange to that extending the parameter space with a point could make a difference as the degenerate case can be approximated arbitrarily well with the non degenerate distribution (except for floating point rounding).
Med venlig hilsen Andreas Noack 2014-09-15 11:25 GMT-04:00 spaceLem <[email protected]>: > Hi John, thanks for your response. > > If you don't think silently accepting lambda = 0 is good practice, how > best should I approach this? The only other thing I can think of is writing > my own function wrapper specifically to catch the case where lambda = 0. > > I don't know how other people work with random numbers, but my suspicions > are that being able to get a random number when lambda = 0 might be more > common than needing skew and kurtosis for that case (although I do accept > these are all problematic cases, ballooning to infinity). The Poisson > distribution is very important to modellers, and zero rates do need to be > handled somehow (currently R, Octave, and the GSL in C++ all accept lambda > = 0, so the textbook definition isn't necessarily the popular one! On the > other hand, Scipy returns an error). > > Regards, > Jamie > > On Monday, 15 September 2014 16:03:32 UTC+1, John Myles White wrote: >> >> I’m not so sure that we should follow the lead of Octave and R here. >> Neither of those languages reify distributions as types, so changes to >> their RNG’s don’t affect other operations on those same distributions. >> >> In contrast, the proposed change here would break a lot of other code in >> Distributions that assumes that every Poisson object defines a non-delta >> distribution. At a minimum, all of the following functions would need to >> have branches inserted into them to work around the lambda = 0 case: >> >> * entropy >> * kurtosis >> * skewness >> >> In addition, every person who ever wrote code in the future that worked >> with Poisson objects would need to know that our definition of the Poisson >> distribution contradicted the definition found in textbooks. This would >> affect people writing code to estimate KL divergences, for example. >> >> — John >> >> On Sep 15, 2014, at 7:33 AM, Andreas Noack <[email protected]> >> wrote: >> >> I can see that R accepts zero as rate parameter so maybe we should do the >> same. Could you open a pull request to Distributions.jl with that change? >> >> Regarding the vectorized version, the answer is that you can do almost >> what you want with a comprehension, i.e. something like X += >> dX*[rand(Poisson(r*dt)) for r in rates]. >> >> Med venlig hilsen >> >> Andreas Noack >> >> 2014-09-15 10:05 GMT-04:00 spaceLem <[email protected]>: >> >>> >>> >>> Hi all, >>> >>> I work in disease modelling, using use a mix of C++ and Octave. I'm >>> fairly new to Julia, although I've managed to get a basic model up and >>> running, and at 1/5.5 times the speed of C++ I'm pretty impressed (although >>> hoping to close in on C++). I'd love to be able to work in one language all >>> the time, and I'm feeling that I'm not far off. >>> >>> I have two questions regarding random numbers and the Poisson >>> distribution. One algorithm I use has a number of possible events, each >>> with an associated event rate. From these events, you choose a time step >>> dt, then the number of times each event happens is Poisson distributed with >>> lambda = rate of the event * dt. In Octave I could write code along these >>> lines (simplified to get the gist of things): >>> >>> rates = 50*rand(1,6); >>> rates(3) = 0; >>> dt = 0.1; >>> K = poissrnd(rates*dt); % = [1 6 0 4 3 4] >>> >>> where K is an array giving the number of times each event occurs. In >>> Julia, I would write >>> >>> using Distributions >>> rates = 50rand(6) >>> rates[3] = 0 >>> dt = 0.1 >>> K = zeros(Float64,6) >>> for i = 1:6; K[i] = rand(Poisson(rates[i] * dt)); end >>> >>> This gives: ERROR: lambda must be positive >>> in Poisson at /Users/spacelem/.julia/Distributions/src/univariate/ >>> poisson.jl:4 >>> in anonymous at no file:1 >>> >>> Which brings me to my first question: how best to handle when the events >>> have zero rates (which is not uncommon, and needs to be handled)? The >>> correct behaviour is that an event with zero rate occurs zero times. >>> >>> I found that editing poisson.jl (mentioned in the error), and changing >>> line 4 from if l > 0 to if l >= 0, then the error goes away and when events >>> have zero rates, they correctly occur zero times. I know the Poisson >>> distribution is technically only defined for lambda > 0, but it really does >>> make sense to handle the case lambda = 0 as returning 0. Somehow I feel >>> that editing that file was probably not the correct thing to do (although >>> it was great that I was able to), but I'd like to follow good practice, and >>> I'm going to run into problems if I ever need to share my code. >>> >>> And my second question: in Octave I can specify an array of lambdas and >>> get back an array of random numbers distributed according to those event >>> rates. Is it possible to do the same in Julia? (I can write rand(6) and get >>> a vector of uniformly distributed random numbers, but I need the for loop >>> to do the same for other distributions). In Octave you can pretty much >>> write something like X += dX * poissrnd(rates * dt); all in one line (where >>> X is a vector of populations, and dX is the event rate / state change >>> matrix), and it would be nice to be as elegant as that in Julia. >>> >>> Thank you very much in advance. >>> >>> Regards, >>> Jamie >>> >> >> >>
