scipy rv_continuous очень медленно

Question

Я использую настраиваемую функцию f(x) для определения настраиваемого распределения с использованием класса copy rv_continuous. Мой код:

class my_pdf_gen(rv_continuous):
    def _pdf(self, x, integral):
        return f(x)/integral

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

my_pdf = my_pdf_gen(my_int,a = a, b = b, name = 'my pdf')

с a,b верхним и нижним пределом диапазона значений и my_int= scipy.integrate.quad(f, a, b)[0]. Я также могу создать случайную выборку данных, используя my_pdf.rvs(my_int, size = 5), но это очень медленно. (До 6 секунд, когда size=9).

Я читал, что следует также перезаписать некоторые другие методы в классе (например, _ppf), но из примеров, которые я обнаружил, мне это не ясно как этого добиться в моем случае.

Большое спасибо!

Answer 1

Я решил проблему, изменив подход и используя метод выборки отбраковки Монте-Карло

def rejection_sampler(p,xbounds,pmax):
    while True:
        x = np.random.rand(1)*(xbounds[1]-xbounds[0])+xbounds[0]
        y = np.random.rand(1)*pmax
        if y<=p(x):
            return x

где p - функция плотности вероятности, xbounds - кортеж, содержащий верхний и нижний пределы pdf, а pmax - максимальное значение pdf в домене.

Здесь был предложен метод отбраковки Монте-Карло: python: случайная выборка из самоопределяемой функции вероятности

Answer 2

Ожидается, что он будет медленным, поскольку реализация generi c выполняет root -решение для cdf, которое само использует численное интегрирование.

Так что лучше всего предоставить _ppf или _rvs реализация. Как это сделать, во многом зависит от деталей f(x). Если вы не можете решить f(x) = r аналитически, рассмотрите возможность табулирования / обратной интерполяции или отклонения выборки.