Наличие N элемента вывода из выборки отклонения для N элементов

Question

Я применяю выборку отклонения для N элементов с учетом функции плотности вероятности pdf. При применении этого метода для N элементов, скорее всего, вы вернете массив значений, который имеет меньшее количество элементов по сравнению с числом N, которое вы оцениваете, то есть от применения метода отклонения без зацикливания для значений condition, которые False.

Чтобы согласовать это, я могу попробовать зациклить значения, которые не соответствуют condition, пока они не станут True. Тем не менее, я не уверен, как зацикливать условие, пока количество элементов в моем массиве не будет равно длине с количеством значений N, учитывая функции, которые я определил.

import numpy

N = 1000 # number of elements
x = np.linspace(0, 200, N) 
pdf = pdf(x) # some pdf


# Rejection Method #1
# -------------------
fx = np.random.random_sample(size=N) * x.max() # uniform random samples scaled out
u = np.random.random_sample(size=N) # uniform random sample
condition = np.where(u <= pdf/pdf.max())[0] # Run first rejection criterion that returns bool values
x_arr = fx[condition] 

# Here, len(x_arr) < N, so I want to fix it until len(x_arr) == N
while len(x_arr) < N:
    ...
    if len(x_arr) == N:
        break

После этого у меня возникают проблемы с формированием метода итерации до len(x_arr) = N.

Answer 1

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

Пример выборки и принятия / отклонения функций:

def sample(N):
    return np.random.uniform(-3, 3, (N,))

def accept(v):
    return np.random.rand(v.size) < stats.norm().pdf(v)

Основной цикл:

def draw(N, f_sample, f_accept):
    out = f_sample(N)
    mask = f_accept(out)
    reject, = np.where(~mask)
    while reject.size > 0:
        fill = f_sample(reject.size)
        mask = f_accept(fill)
        out[reject[mask]] = fill[mask]
        reject = reject[~mask]
    return out

Проверка работоспособности:

>>> counts, bins = np.histogram(draw(100000, sample, accept))
>>> counts / stats.norm().pdf((bins[:-1] + bins[1:]) / 2)
array([65075.50020815, 65317.17811578, 60973.84255365, 59440.53739031,
       58969.62310004, 59267.33983256, 60565.1928325 , 61108.60840388,
       64303.2863583 , 68293.86441234])

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