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 >
