As this situation arises in other cases, I've opened an issue here: https://github.com/JuliaStats/Distributions.jl/issues/283
The skewness shouldn't be Inf as it arises as 0/0. We could return a NaN though. On Tuesday, 16 September 2014 12:30:20 UTC+1, spaceLem wrote: > > I suppose it comes from one's perspective -- from a modeller's point of > view, a zero rate is certainly not a special case! The error was > surprising, which made me baulk at handling it. > > From a purity argument, I agree with you -- Poisson(0) is not defined, > however from a usefulness argument, simply returning 0 is what people > (especially those coming from most other languages) would expect (after > all, the limit of P(0) as lambda -> 0 is 0), and returning infinity for > skewness etc. rather than an error also makes a geometric sense. > > However, I'll probably handle it as David Gonzales' nice one liner for the > moment. > > Thanks for your advice! > > Regards, > Jamie > > On Monday, 15 September 2014 16:43:13 UTC+1, John Myles White wrote: >> >> Why not just write out an explicit check in your code for the special >> case? Here's an example of how one might do that: >> >> n_samples = 10 >> n_lambdas = 3 >> lambdas = [0.0, 0.1, 0.2] >> >> samples = Array(Int, n_samples, n_lambdas) >> >> for (i, lambda) in enumerate(lambdas) >> for sample in 1:n_samples >> if rate == 0.0 >> samples[sample, i] = 0 >> else >> samples[sample, i] = rand(Poisson(lambda)) >> end >> end >> end >> >> In my mind, this is the best possible approach to this kind of problem: >> you want slighly unusual behavior, so you should make the place where your >> goals are slightly unusual as explicit as possible so that anyone who reads >> your code will understand exactly what you're doing and exactly where >> you're breaking from textbook definitions. >> >> As for precedent, I believe the list of systems you've provided breaks >> down as follows: >> >> RNG's take in numbers, not objects: >> >> * R >> * Octave >> * GSL >> >> RNG's take in objects, not numbers: >> >> * SciPy >> * Julia >> >> So it seems like the split in behavior is perfectly explicable in terms >> of the types of arguments you provide to RNG's. This seems like a very good >> argument for not changing the behavior of Distributions. >> >> -- John >> >> On Sep 15, 2014, at 8:25 AM, spaceLem <[email protected]> wrote: >> >> 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 >>>> >>> >>
