scipy.stats.sampling.

RatioUniforms#

класс scipy.stats.sampling.RatioUniforms(pdf, *, umax, vmin, vmax, c=0, random_state=None)[источник]#

Генерация случайных выборок из функции плотности вероятности с использованием метода отношения равномерных распределений.

Параметры:

pdfcallable: Функция с сигнатурой pdf(x) которая пропорциональна функции плотности вероятности распределения.
umaxfloat: Верхняя граница ограничивающего прямоугольника в u-направлении.
vminfloat: Нижняя граница ограничивающего прямоугольника в направлении v.
vmaxfloat: Верхняя граница ограничивающего прямоугольника в v-направлении.
cfloat, опционально.: Параметр сдвига метода отношения равномерных величин, см. Примечания. По умолчанию 0.
random_state{None, int, numpy.random.Generator,: numpy.random.RandomState, опционально

Если seed равно None (или np.random), numpy.random.RandomState используется синглтон. Если seed является int, новый RandomState используется экземпляр, инициализированный с seed. Если seed уже является Generator или RandomState экземпляр, тогда этот экземпляр используется.

Методы

rvs([size])

Выборка случайных величин

Примечания

Дана одномерная функция плотности вероятности pdf и константа c, определяют множество A = {(u, v) : 0 < u <= sqrt(pdf(v/u + c))}. Если (U, V) является случайным вектором, равномерно распределенным по A, затем V/U + c следует распределению согласно pdf.

Вышеуказанный результат (см. [1], [2]) можно использовать для выборки случайных величин используя только PDF, т.е. инверсия CDF не требуется. Типичные выборы c равны нулю или мода pdf. Множество A является подмножеством прямоугольника R = [0, umax] x [vmin, vmax] где

umax = sup sqrt(pdf(x))
vmin = inf (x - c) sqrt(pdf(x))
vmax = sup (x - c) sqrt(pdf(x))

В частности, эти значения конечны, если pdf ограничена и x**2 * pdf(x) ограничена (т.е. субквадратичные хвосты). Можно сгенерировать (U, V) равномерно на R и возвращает V/U + c if (U, V) также находятся в A который можно напрямую проверить.

Алгоритм не изменяется, если заменить pdf на k * pdf для любой константы k > 0. Таким образом, часто удобно работать с функцией, пропорциональной функции плотности вероятности, отбрасывая ненужные нормировочные множители.

Интуитивно, метод работает хорошо, если A заполняет большую часть ограничивающего прямоугольника так, что вероятность высока, что (U, V) лежит в A всякий раз, когда он лежит в R поскольку количество требуемых итераций становится слишком большим в противном случае. Чтобы быть более точным, заметим, что ожидаемое количество итераций для выборки (U, V) равномерно распределены на R такой, что (U, V) также находится в A задается отношением area(R) / area(A) = 2 * umax * (vmax - vmin) / area(pdf), где area(pdf) является интегралом от pdf (что равно единице, если используется функция плотности вероятности, но может принимать другие значения, если используется функция, пропорциональная плотности). Равенство выполняется, поскольку площадь A равно 0.5 * area(pdf) (Теорема 7.1 в [1]). Если выборка не может сгенерировать ни одной случайной величины после 50000 итераций (т.е. ни один выбор не находится в A), возникает исключение.

Если ограничивающий прямоугольник указан некорректно (т.е. если он не содержит A), алгоритм выбирает из распределения, отличного от заданного pdf. Поэтому рекомендуется выполнить тест, такой как kstest в качестве проверки.

Ссылки

[1] (1,2)

L. Devroye, «Non-Uniform Random Variate Generation», Springer-Verlag, 1986.

[2]

W. Hoermann и J. Leydold, «Generating generalized inverse Gaussian random variates», Statistics and Computing, 24(4), стр. 547–557, 2014.

[3]

A.J. Kinderman и J.F. Monahan, "Генерация случайных величин на компьютере с использованием отношения равномерных отклонений", ACM Transactions on Mathematical Software, 3(3), стр. 257–260, 1977.

Примеры

>>> import numpy as np
>>> from scipy import stats

>>> from scipy.stats.sampling import RatioUniforms
>>> rng = np.random.default_rng()

Моделировать нормально распределённые случайные величины. В этом случае легко вычислить ограничивающий прямоугольник явно. Для простоты мы опускаем нормировочный множитель плотности.

>>> f = lambda x: np.exp(-x**2 / 2)
>>> v = np.sqrt(f(np.sqrt(2))) * np.sqrt(2)
>>> umax = np.sqrt(f(0))
>>> gen = RatioUniforms(f, umax=umax, vmin=-v, vmax=v, random_state=rng)
>>> r = gen.rvs(size=2500)

K-S тест подтверждает, что случайные величины действительно нормально распределены (нормальность не отвергается на 5% уровне значимости):

>>> stats.kstest(r, 'norm')[1]
0.250634764150542

Экспоненциальное распределение предоставляет ещё один пример, где ограничивающий прямоугольник может быть определён явно.

>>> gen = RatioUniforms(lambda x: np.exp(-x), umax=1, vmin=0,
...                     vmax=2*np.exp(-1), random_state=rng)
>>> r = gen.rvs(1000)
>>> stats.kstest(r, 'expon')[1]
0.21121052054580314