Thanks for your response! I can do that (once I figure out how). Regards, Jamie
On Monday, 15 September 2014 15:33:25 UTC+1, Andreas Noack 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 <spac...@minioak.com <javascript:>>: > >> >> >> 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 >> > >
