Я пытаюсь написать программу на R, которая имитирует псевдослучайные числа из распределения с помощью кумулятивной функции распределения:
где
Я пробовал выборку с обратным преобразованием, но обратное не представляется аналитически разрешимым. Я был бы рад, если бы вы могли предложить решение этой проблемы
r
random-generation
Себастьян
источник
источник
Ответы:
Существует простое (и, если я могу добавить, элегантное) решение этого упражнения: поскольку1−F(x) выглядит как произведение двух распределений выживания:
(1−F(x))=exp{−ax−bp+1xp+1}=exp{−ax}1−F1(x)exp{−bp+1xp+1}1 - F2( х )
распределение F - это распределение
Икс= min { X1, X2}Икс1∼ F1, X2∼ F2
В этом случае F1 является экспоненциальным распределением Е( а ) а F2 является 1 / ( р + 1 ) -ым степень экспоненциального распределения Е( б / ( р + 1 ) ) .
Соответствующий код R настолько прост, насколько это возможно
и это определенно намного быстрее, чем обратное разрешение pdf и accept-reject:
с удивительно идеальной посадкой:
источник
Вы всегда можете численно решить обратное преобразование.
Ниже я делаю очень простой поиск пополам. Для заданной входной вероятности (я использую так как в вашей формуле уже есть ), я начну с и . Затем я удваиваю до пор, пока . Наконец, я итеративно пополам интервал до тех пор, пока его длина не станет меньше а его средняя точка удовлетворит .Q Q п ИксL= 0 Икср= 1 Икср F( хр) > д [ хL, хр] ε ИксM F( хM) ≈ q
ECDF подходит вашему достаточно хорошо для моих вариантов и , и это достаточно быстро. Вероятно, вы могли бы ускорить это, используя некоторую оптимизацию типа Ньютона вместо простого поиска по разделению пополам.F a b
источник
Существует несколько запутанная, если прямое решение путем принятия-отклонения. Во-первых, простое дифференцирование показывает, что pdf распределения - это Во-вторых, так как у нас есть верхняя граница третьих, учитывая второе слагаемое в , возьмем замену переменной , т. е. . Тогда является якобианом изменения переменной. Еслие( х ) = ( а + б хп) опыт{ - а х - бр + 1Икср + 1} е( х ) = а е- хе- б хр + 1/ (р+1)≤ 1+ б хпе- б хр + 1/ (р+1)е- х≤ 1 е( х ) ≤ г( х ) = а е- х+ б хпе- б хр + 1/ (р+1) грамм ξ= хр + 1 х = ξ1 / ( р + 1 ) д хд ξ= 1р + 1ξ1р + 1- 1= 1р + 1ξ- рр + 1 Икс имеет плотность вида где - нормализующая константа, тогда имеет плотность
что означает , что (я) является распределяется как экспоненциальная переменная и (ii) константа равна единице. Следовательно, оказывается равным взвешенной смеси экспоненциального распределения и -й степени экспоненциальногоκ b xпе- б хр + 1/ (р+1) κ Ξ = X1 / ( р + 1 ) κ b ξпр + 1е- b ξ/ (р+1)1р + 1ξ- рр + 1= κ бр + 1е- b ξ/ (р+1) Ξ Е(b/(p+1)) κ g(x) E( а ) 1 / ( р + 1 ) Е( б/ ( р + 1 ) ) распределение, по модулю отсутствующей мультипликативной константы для учета весов:
И легко смоделировать как смесь.2 е(х ) ≤ г( х ) = 2 ( 12е- ах+ 12б хпе- бхр + 1/ (р+1)) грамм
Таким образом, R-рендеринг алгоритма принятия-отклонения
и для n-образца:
Вот иллюстрация для a = 1, b = 2, p = 3:
источник