Implement dist_rand for dislnorm and disexp

[linzhp]

@csgillespie

Please review. Specifically, please confirm my use of m$pars in disexp$dis_rand is correct.

Thanks!

[csgillespie]

Thanks for the request.

A couple of suggestions/comments before merging:

• Your use of m$pars is correct
• Don't grow your vector as this is very inefficient. So have result = numeric(n)
• In the lognormal RNG, it should be log(m$xmin -0.5) since we are discretising the distribution.
• For the exponential random numbers, a non-rejection algorithm is more efficient. I think,

u = runif(n, 0, exp(-m$pars*(m$xmin-0.5)))
-log(u)/m*pars

will work - I just rearranged the standard exponential random number algorithm.

What do you think?

[linzhp]

Thanks for the comments. I will address them and update the pull request.

Question: in the lognormal RNG, wouldn't min <- round(log(m$xmin)) be better than min <- log(m$xmin -0.5)?

[csgillespie]

I think you're correct about rounding xmin. A couple of other points I just thought of:

1. Check your RNGs by simulating a large number of times and comparing to the exact pdf - dist_pdf
2. For the lognormal example, if we simulate n values on average we will discard plnorm(xmin-0.5, m$pars[1],m$pars[2], lower.tail=TRUE) So simulate (something like)

ceiling(n/plnorm(m$xmin-0.5,  m$pars[1],m$pars[2], lower.tail=FALSE))

[linzhp]

I've updated the pull request according to your suggestions. However, I don't understand what

ceiling(n/plnorm(m$xmin-0.5,  m$pars[1],m$pars[2], lower.tail=FALSE))

does.

I used

r<-dist_rand(exponential, 1000)
plot(density(r, from=min(r)))
curve(dexp(x, exponential$pars), from=min(r), to=max(r))

to visually compare the two pdfs