Как генерировать случайные числа из отсечения log-normal распределения в 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

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

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

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

Answer 2

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

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

Answer 3

Возможно, излишний для вашего распространения - но вы всегда можете написать Metropolis sampler.

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

Answer 4

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

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

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

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

y = rand (size_dist,1); 

с последующим способом, чтобы поместить распределение вдоль КОР. Любая техника будет работать, но здесь это самый простой (хотя и безвкусный)

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)