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] <javascript:>> > 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 >>> >> >
