How to generate random numbers from cut-off log-normal distribution in Matlab?

Question

The radii r is drawn from a cut-off log-normal distribution, which has a following probability density function:

pdf=((sqrt(2).*exp(-0.5*((log(r/rch)).^2)))./((sqrt(pi.*(sigma_nd.^2))...
    .*r).*(erf((log(rmax/rch))./sqrt(2.*(sigma_nd.^2)))-erf((log(rmin/rch))./sqrt(2.*(sigma_nd.^2))))));

rch, sigma_nd, rmax, and rmin are all constants.

I found the explanation from the net, but it seems difficult to find its integral and then take inverse in Matlab.

Answer 1

I checked, but my first instinct is that it looks like log(r/rch) is a truncated normal distribution with limits of log(rmin/rch) and log(rmax/rch). So you can generate the appropriate truncated normal random variable, say y, then r = rch * exp(y).

You can generate truncated normal random variables by generating the untruncated values and replacing those that are outside the limits. Alternatively, you can do it using the CDF, as described by @PengOne, which you can find on the wikipedia page.

I'm (still) not sure that your p.d.f. is completely correct, but what's most important here is the distribution.

Answer 2

If your PDF is continuous, then you can integrate to get a CDF, then find the inverse of the CDF and evaluate that at the random value.

If your PDF is not continuous, then you can get a discrete CDF using cumsum, and use that as your initial Y value in interp(), with the initial X value the same as the values the PDF was sampled at, and asking to interpolate at your array of rand() numbers.

Answer 3

it's not clear what's your variable, but i'm assuming it's r.

the simplest way to do this is, as Chris noted, first get the cdf (note that if r starts at 0, pdf(1) is Nan... change it to 0):

cdf = cumtrapz(pdf);
cdf = cdf / cdf(end);

then spawn a uniform distribution (size_dist indicating the number of elements):

y = rand (size_dist,1);

followed by a method to place distribution along the cdf. Any technique will work, but here is the simplest (albeit inelegant)

x = zeros(size_dist,1);
for i = 1:size_dist
    x(i) = find( y(i)<= cdf,1);
end

and finally, returning to the original pdf. Use matlab numerical indexing to reverse course. Note: use r and not pdf:

pdfHist = r(x);
hist (pdfHist,50)

Answer 4

Probably an overkill for your distribution - but you can always write a Metropolis sampler.

On the other hand - implementation is straight forward so you'd have your sampler very quick.