Как правильно пробовать плотность?

Question

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

Вот пример, где я пытался образец нормального распределения:

n = 100000 
xx = np.random.uniform(-5, 5, n) 
rho = mpl.pylab.normpdf(xx, 0, 1) 
rnd = np.random.rand(n) 
ix = np.where(rho > rnd) 
xx = xx[ix] 
h = plt.hist(xx, bins=20, normed=True) 
# plot density 
x = np.linspace(-5, 5, 100) 
plt.plot(x, mpl.pylab.normpdf(x, 0, 1)) 

Он работает и Я получил:

Теперь, если я изменил плотность, я не правильно пробовать его. Я проверил, хорошо ли нормируется плотность, и это так. Таким образом, я не понимаю, где я не прав

n = 100000 
xx = np.random.uniform(0, 1, n) 
rho = 2 * np.sin(2 * xx * np.pi)**2 
rnd = np.random.rand(n) 
ix = np.where(rho > rnd) 
xx = xx[ix] 
h = plt.hist(xx, bins=20, normed=True) 
# plot density 
x = np.linspace(0, 1, 100) 
print(np.trapz(2 * np.sin(2 * x * np.pi)**2, x)) 
plt.plot(x, 2 * np.sin(2 * x * np.pi)**2) 

Answer 1

Вы делаете rejection sampling

В первом случае максимальное значение PDF является < 1, и вы рисуете rnd из [0,1] , поэтому все значения ниже макс. Вы выбрасывая больше значения, чем нужно, хотя, так как максимум строго меньше 1. Во втором случае максимальная ПРВ является но вы все еще рисунок rnd из [0,1] в линии

rnd = np.random.rand(n) 

Вы должны изменить эту строку, чтобы она равномерно от [0,2]. Обратите внимание, что несколько плоских вершин ваших гистограмм соответствуют частям [0,1], где pdf> 1. У вашего кода нет способа рассматривать некоторые из этих значений иначе, чем другие.

Answer 2

В первом примере вы отвергаете слишком много, а во втором - недостаточно.

Оптимальный чехол, когда вы отбираете Y от 0 до PDF _max.

В первом случае будем называть

rnd = np.random.rand(n)/np.sqrt(2.0 * np.pi) 

Во втором случае

rnd = 2.0 * np.random.rand(n)