генерация случайной выборки из экспоненциального распределения в Stata

Я пытаюсь выполнить моделирование в Stata со случайной выборкой из 10000 для (i) переменной X с pdf f(x) = 2*x*exp(-x^2), X>0 и (ii) Y=X^2. Я вычислил cdf F как 1-exp(-x^2), поэтому инверсия F равна sqrt(-ln(1-u). I использовал следующий код в Stata:

(1)  
 clear  
 set obs 10000  
 set seed 527665  
 gen u= runiform()  
 gen x= sqrt(-ln(1-u))  
 histogram x  
 summ x, detail  
(mean 0.88, sd 0.46)  
  

(2)  
clear  
set obs 10000  
set seed 527665  
gen u= runiform()  
gen x= (sqrt(-ln(1-u)))^2  
summ x, detail  
(mean 0.99, sd 0.99) 

(3)    
clear  
set obs 10000  
set seed 527665  
gen u= rexponential(1)  
gen x= 2*u*exp(-(u^2))  
summ x, detail  
(mean 0.49, sd 0.28)  

(4)
clear  
set obs 10000  
set seed 527665  
gen v= runiform()  
gen u=1/v  
gen x= 2*u*exp(-(u^2))  
histogram x  
summ x, detail  
(mean 0.22, sd 0.26)

Мои запросы: (i) (1) и (2) основаны на преобразовании интеграла вероятности, с которым я столкнулся, но не понимаю. Если (1) и (2) являются допустимыми подходами, то какая интуиция стоит за этим, (ii) вывод для (3) не кажется правильным; Я не уверен, правильно ли я применяю реэкспоненциальную функцию и каков параметр масштаба (похоже, в справке по этому поводу нет объяснения) (iii) вывод для (4) также не кажется правильным, и я был интересно, почему этот подход ошибочен.

Спасибо


random stata exponential-distribution
person Matt    schedule 11.09.2020    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Что ж, то, что вы придумали в качестве дистрибутива, мне кажется нормальным.

If

PDF(x) = 2 x exp(-x2), x в [0...бесконечность), тогда

CDF(x) = 1 - exp(-x2)

что означает, что это в основном квадратный корень из экспоненциально распределенного RV. Экспоненциальное распределение выборка выполняется с использованием -ln(1-u) или -ln(u)

У меня нет Stata, просто смотрю код

(1) выглядит нормально, вы сэмплируете экспоненту и получаете из нее квадратный корень

(2) похоже, что вы выбираете квадратный корень из экспоненты и сразу же возводите его обратно. Вы вернетесь экспоненциально, я верю

(3) Я не знаю, что это должно означать, показатель степени квадрата экспоненты? Должно быть

clear  
set obs 10000  
set seed 527665  
gen u = rexponential(1)  
gen x = sqrt(u)
summ x, detail

reexponential() — это то же самое, что и -ln(1-runiform())

(4) Не имеет смысла. Показатель от квадрата униформы?

Я быстро написал простой код Python для иллюстрации

import numpy as np
import matplotlib.pyplot as plt

x = np.random.random(100000) // uniform in [0...1)
xx = np.sqrt(-np.log(1.0-x)) // -log(1-x) is exponential, then square root

q = np.linspace(0.0, 3.0, 101)
z = 2.0*q*np.exp(-q*q)

n, bins, patches = plt.hist(xx, 50, density=True, facecolor='g', alpha=0.75)
plt.plot(q, z, 'r-')
plt.show()

с изображением

введите здесь описание изображения

person Severin Pappadeux    schedule 11.09.2020
comment
Для (4), если gen v= runiform() генерирует случайные числа от 0 до 1, будет ли функция, обратная v, 1/v, давать случайные числа от 0 до бесконечности - person Matt; 11.09.2020
comment
@Matt хорошо, и это все еще не имеет смысла. В вашем PDF нет разделения, почему вы хотите разделить равномерно распределенные RV? Это должен быть квадратный корень из экспоненты либо через -ln(1-u), либо через rexponential(), как я указал в коде. Вы можете нарисовать гистограммы и сравнить их - person Severin Pappadeux; 11.09.2020
comment
@Matt Я обновил ответ, добавил перекрывающийся PDF поверх гистограммы, все выглядит хорошо - person Severin Pappadeux; 11.09.2020
comment
gen x = sqrt(rexponential(1)) немного более прямой. - person Nick Cox; 11.09.2020
comment
@NickCox конечно, наверное, немного эффективнее - person Severin Pappadeux; 11.09.2020