Would be good to sort out the implications, but I'm hoping that atomic measures are not as problematic as JMW fears.
On Tuesday, September 16, 2014 8:14:58 AM UTC-4, Simon Byrne wrote: > > 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 >>>>> >>>> >>>
