Как сгенерировать случайные числа из логарифмического нормального распределения в Matlab?

Question

Радиусы r взяты из нормального логарифмического распределения, которое имеет следующую функцию плотности вероятности:

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 и rmin - все константы.

Я нашел объяснение в сети, но, кажется, трудно найти его интеграл, а затем взять обратное в Matlab.

Answer 1

не ясно, какая у вас переменная, но я предполагаю, что это r.

Самый простой способ сделать это, как заметил Крис, сначала получить cdf (обратите внимание, что если r начинается с 0, pdf(1) равно Nan ... измените его на 0):

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

затем создает равномерное распределение (size_dist, указывающее количество элементов):

y = rand (size_dist,1);

, за которым следует метод для распространения по cdf. Любая техника будет работать, но здесь самый простой (хотя и не элегантный)

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

и, наконец, возвращаемся к исходному pdf. Используйте числовое индексирование Matlab, чтобы изменить курс. Примечание: используйте r, а не pdf:

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

Answer 2

Я проверил, но мой первый инстинкт заключается в том, что log(r/rch) - это усеченное нормальное распределение с пределами log(rmin/rch) и log(rmax/rch).Таким образом, вы можете сгенерировать соответствующую усеченную нормальную случайную переменную, скажем y, затем r = rch * exp(y).

. Вы можете сгенерировать усеченные нормальные случайные переменные, сгенерировав неусеченные значения и заменив их значения, выходящие за пределы.Кроме того, вы можете сделать это с помощью CDF, как описано в @PengOne, которое вы можете найти на странице википедии .

Я (все еще) не уверен, что ваш pdf полностьюправильно, но самое главное здесь - это дистрибутив.

Answer 3

Если ваш PDF-файл является непрерывным, то вы можете интегрировать его, чтобы получить CDF, а затем найти инверсию CDF и оценить ее по случайному значению.

Если ваш PDF-файл не является непрерывным, тогда вы можете получитьдискретный CDF, использующий cumsum, и используйте его в качестве начального значения Y в interp (), причем начальное значение X совпадает со значениями, в которых был выбран файл PDF, и запрашивает интерполяцию в массиве чисел rand ().

Answer 4

Возможно, это избыточное решение для вашего дистрибутива, но вы всегда можете написать Метрополис сэмплер .

С другой стороны - реализация проста, поэтому ваш сэмплер будет очень быстрым.