or with comprehension and unicode-sugar: using Distributions Λ = 50rand(6) ; Λ[3] = 0.0 ; Δt = 0.1 [λ > 0.0 ? rand(Poisson(λ*Δt)) : 0 for λ in Λ]
@assert Julia > max(Matlab,Octave,c++) On Monday, September 15, 2014 6:51:39 PM UTC+3, John Myles White wrote: > > The KL divergence (and other comparisons) between this edge case and any > standard instance of the Poisson distribution are problematic because the > inclusion of this point in parameter space changes the support of the > distribution from the entire set of non-negative integers to the set > containing only zero. > > -- John > > On Sep 15, 2014, at 8:42 AM, Andreas Noack <[email protected] > <javascript:>> wrote: > > This is not a big deal to me, but what exactly would break by allowing > λ=0? Returning inf for skewness and kurtosis and an error for entropy > doesn't seem problematic to me. I don't do estimation of KL divergences, > but it appears a bit strange to that extending the parameter space with a > point could make a difference as the degenerate case can be approximated > arbitrarily well with the non degenerate distribution (except for floating > point rounding). > > Med venlig hilsen > > Andreas Noack > > 2014-09-15 11:25 GMT-04:00 spaceLem <[email protected] <javascript:>>: > >> 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 >>>> >>> >>> >>> > >
