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
