Полиномиальная интерполяция на основе инверсии CDF (PINV)

• Требуется: PDF

• Опционально: CDF, мода, центр

• Скорость:

  • Настройка: (очень) медленная

  • Сэмплирование: (очень) быстрое

Полиномиальная интерполяция на основе инверсии CDF (PINV) — это метод инверсии, который требует только функцию плотности для выборки из распределения. Он основан на полиномиальной интерполяции PPF и интегрировании Гаусса-Лобатто PDF. Он обеспечивает контроль над ошибкой интерполяции и ошибкой интегрирования. Его основная цель — обеспечить очень быструю выборку, которая почти одинакова для любого заданного распределения за счёт умеренного или медленного времени настройки. Это самый быстрый известный метод инверсии для случая с фиксированными параметрами.

Метод инверсии является самым простым и гибким для выборки неоднородных случайных величин. Для целевого распределения с функцией распределения \(F\) и равномерная случайная величина \(U\) выбранный из \(\text{Uniform}(0, 1)\), случайная вариация X генерируется путем преобразования равномерной случайной вариации \(U\) используя PPF (обратную функцию распределения) распределения:

\[X = F^{-1}(U)\]

Этот метод подходит для стохастического моделирования благодаря своим преимуществам. Некоторые из наиболее привлекательных:

• Он сохраняет структурные свойства генератора равномерных случайных чисел.

• Преобразует равномерную случайную величину \(U\) взаимно однозначно в неравномерные случайные величины \(X\).

• Простой и эффективный выборка из усеченных распределений.

К сожалению, PPF редко доступен в замкнутой форме или слишком медленный, когда доступен. Для многих распределений CDF также нелегко получить. Этот метод решает обе проблемы. Пользователю нужно только предоставить PDF и, опционально, точку около моды (называемую «центром») вместе с размером максимально допустимой ошибки. Он использует комбинацию адаптивной и простой квадратуры Гаусса-Лобатто для получения CDF и интерполяции Ньютона для получения PPF. Метод не является точным, так как он производит только случайные величины аппроксимированного распределения. Тем не менее, максимально допустимая ошибка аппроксимации может быть установлена близко к машинной точности. Для измерения и контроля ошибки используется концепция u-ошибки, которая определяется как:

\[\epsilon_{u}(u) = | u - F\left(F^{-1}_{a}(u)\right) |\]

где \(u \in (0, 1)\) — это квантиль, где мы хотим измерить ошибку и \(F^{-1}_a\) является аппроксимированной PPF данного распределения.

Максимальная u-ошибка — это критерий для ошибок аппроксимации при численном вычислении CDF и PPF. Максимально допустимая u-ошибка алгоритма называется u-разрешением алгоритма и обозначается \(\epsilon_{u}\):

\[\sup_{u \in (0,1)} | u - F\left(F^{-1}_{a}(u)\right) | \le \epsilon_{u}\]

Основное преимущество u-ошибки в том, что её можно легко вычислить, если доступна CDF. Мы ссылаемся на [1] для более подробного обсуждения.

Кроме того, метод работает только для ограниченных распределений. В случае бесконечных хвостов, концы хвостов обрезаются так, что площадь под ними меньше или равна \(0.05\epsilon_{u}\).

Существуют некоторые ограничения для данного распределения:

• Носитель распределения (т.е. область, где плотность вероятности строго положительна) должен быть связным. На практике это означает, что область, где плотность вероятности «не слишком мала», должна быть связной. Унимодальные плотности удовлетворяют этому условию. Если это условие нарушено, то область определения распределения может быть усечена.

• При численном интегрировании PDF, заданная функция плотности вероятности должна быть непрерывной и гладкой.

• PDF должен быть ограниченным.

• Алгоритм имеет проблемы, когда распределение имеет тяжелые хвосты (тогда обратная CDF становится очень крутой при 0 или 1) и требуемое разрешение u очень мало. Например, распределение Коши, вероятно, покажет эту проблему, когда требуемое разрешение u меньше 1.e-12.

Следующие четыре шага выполняются алгоритмом во время настройки:

• Вычисление конечных точек распределения: если задана конечная поддержка, этот шаг пропускается. В противном случае концы хвостов обрезаются так, чтобы площадь под ними была меньше или равна \(0.05\epsilon_{u}\).

• Область определения разделена на подинтервалы для вычисления CDF и PPF.

• CDF вычисляется с использованием квадратуры Гаусса-Лобатто, так что ошибка интегрирования не превышает \(0.05I_{0}\epsilon_{u}\) где \(I_{0}\) приблизительно равна общей площади под PDF.

• PPF вычисляется с использованием интерполяционной формулы Ньютона с максимальной ошибкой интерполяции \(0.9\epsilon_{u}\).

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

>>> import numpy as np
>>> from scipy.stats.sampling import NumericalInversePolynomial
>>> class StandardNormal:
...     def pdf(self, x):
...         return np.exp(-0.5 * x*x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> rng = NumericalInversePolynomial(dist, random_state=urng)

Генератор настроен, и мы можем начать выборку из нашего распределения:

>>> rng.rvs((5, 3))
array([[-1.52449963,  1.31933688,  2.05884468],
       [ 0.48883931,  0.15207903, -0.02150773],
       [ 1.11486463,  1.95449597, -0.30724928],
       [ 0.9854643 ,  0.29867424,  0.7560304 ],
       [-0.61776203,  0.16033378, -1.00933003]])

Мы можем посмотреть на гистограмму случайных величин, чтобы проверить, насколько хорошо они соответствуют нашему распределению:

>>> import matplotlib.pyplot as plt
>>> from scipy.stats import norm
>>> from scipy.stats.sampling import NumericalInversePolynomial
>>> class StandardNormal:
...     def pdf(self, x):
...         return np.exp(-0.5 * x*x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> rng = NumericalInversePolynomial(dist, random_state=urng)
>>> rvs = rng.rvs(10000)
>>> x = np.linspace(rvs.min()-0.1, rvs.max()+0.1, num=10000)
>>> fx = norm.pdf(x)
>>> plt.plot(x, fx, "r-", label="pdf")
>>> plt.hist(rvs, bins=50, density=True, alpha=0.8, label="rvs")
>>> plt.xlabel("x")
>>> plt.ylabel("PDF(x)")
>>> plt.title("Samples drawn using PINV method.")
>>> plt.legend()
>>> plt.show()

Максимально допустимая ошибка (т.е. u-разрешение) может быть изменена путем передачи u_resolution ключевое слово при инициализации:

>>> rng = NumericalInversePolynomial(dist, u_resolution=1e-12,
...                                  random_state=urng)

Это приводит к более точному приближению PPF, и сгенерированные СВ более точно следуют точному распределению. Однако обратите внимание, что это происходит за счет дорогой настройки.

Время настройки в основном зависит от количества вычислений PDF. Оно выше для PDF, которые сложно вычислять. Заметим, что мы можем игнорировать константу нормализации для ускорения вычислений PDF. Вычисления PDF увеличиваются во время настройки для малых значений u_resolution.

>>> from scipy.stats.sampling import NumericalInversePolynomial
>>> class StandardNormal:
...     def __init__(self):
...         self.callbacks = 0
...     def pdf(self, x):
...         self.callbacks += 1
...         return np.exp(-0.5 * x*x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> # u_resolution = 10^-8
>>> # => less PDF evaluations required
>>> # => faster setup
>>> rng = NumericalInversePolynomial(dist, u_resolution=1e-8,
...                                  random_state=urng)
>>> dist.callbacks
4095        # may vary
>>> dist.callbacks = 0  # reset the number of callbacks
>>> # u_resolution = 10^-10 (default)
>>> # => more PDF evaluations required
>>> # => slow setup
>>> rng = NumericalInversePolynomial(dist, u_resolution=1e-10,
...                                  random_state=urng)
>>> dist.callbacks
11454       # may vary
>>> dist.callbacks = 0  # reset the number of callbacks
>>> # u_resolution = 10^-12
>>> # => lots of PDF evaluations required
>>> # => very slow setup
>>> rng = NumericalInversePolynomial(dist, u_resolution=1e-12,
...                                  random_state=urng)
13902     # may vary

Как мы видим, количество требуемых вычислений PDF очень велико, и быстрый PDF критически важен для алгоритма. Однако это помогает уменьшить количество подынтервалов, необходимых для достижения цели по ошибке, что экономит память и делает выборку быстрой. NumericalInverseHermite является аналогичным методом инверсии, который инвертирует CDF на основе интерполяции Эрмита и обеспечивает контроль над максимально допустимой ошибкой через u-разрешение. Но он использует гораздо больше интервалов по сравнению с NumericalInversePolynomial:

>>> from scipy.stats.sampling import NumericalInverseHermite
>>> # NumericalInverseHermite accepts a `tol` parameter to set the
>>> # u-resolution of the generator.
>>> rng_hermite = NumericalInverseHermite(norm(), tol=1e-12)
>>> rng_hermite.intervals
3000
>>> rng_poly = NumericalInversePolynomial(norm(), u_resolution=1e-12)
>>> rng_poly.intervals
252

Когда точная CDF распределения доступна, можно оценить u-ошибку, достигнутую алгоритмом, вызвав u_error method:

>>> from scipy.special import ndtr
>>> class StandardNormal:
...     def pdf(self, x):
...         return np.exp(-0.5 * x*x)
...     def cdf(self, x):
...         return ndtr(x)
...
>>> dist = StandardNormal()
>>> urng = np.random.default_rng()
>>> rng = NumericalInversePolynomial(dist, random_state=urng)
>>> rng.u_error(sample_size=100_000)
UError(max_error=8.785949745515609e-11, mean_absolute_error=2.9307548109436816e-11)

u_error запускает монте-карло симуляцию с заданным количеством выборок для оценки u-ошибки. В приведённом выше примере симуляция использует 100 000 выборок для аппроксимации u-ошибки. Она возвращает максимальную u-ошибку (max_error) и средняя абсолютная u-ошибка (mean_absolute_error) в UError namedtuple. Как мы видим, max_error ниже значения по умолчанию u_resolution (1e-10).

Также возможно вычислить PPF заданного распределения после инициализации генератора:

>>> rng.ppf(0.975)
1.959963985701268
>>> norm.ppf(0.975)
1.959963984540054

Мы можем использовать это, например, для проверки максимальной и средней абсолютной u-ошибки:

>>> u = np.linspace(0.001, 0.999, num=1_000_000)
>>> u_errors = np.abs(u - dist.cdf(rng.ppf(u)))
>>> u_errors.max()
8.78600525666684e-11
>>> u_errors.mean()
2.9321444940323206e-11

Приближенный метод PPF, предоставляемый генератором, вычисляется намного быстрее, чем точный PPF распределения.

Во время настройки сохраняется таблица точек CDF, которая может использоваться для аппроксимации CDF после создания генератора:

>>> rng.cdf(1.959963984540054)
0.9750000000042454
>>> norm.cdf(1.959963984540054)
0.975

Мы можем использовать это, чтобы проверить, превышает ли ошибка интегрирования при вычислении CDF \(0.05I_{0}\epsilon_{u}\). Здесь \(I_0\) является \(\sqrt{2\pi}\) (нормировочная константа для стандартного нормального распределения):

>>> x = np.linspace(-10, 10, num=100_000)
>>> x_error = np.abs(dist.cdf(x) - rng.cdf(x))
>>> x_error.max()
4.506062190046123e-12
>>> I0 = np.sqrt(2*np.pi)
>>> max_integration_error = 0.05 * I0 * 1e-10
>>> x_error.max() <= max_integration_error
True

Таблица CDF, вычисленная во время настройки, используется для оценки CDF, и требуется лишь некоторая дополнительная тонкая настройка. Это уменьшает вызовы PDF, но поскольку шаг тонкой настройки использует простую квадратуру Гаусса-Лобатто, PDF вызывается несколько раз, замедляя вычисления.

Ссылки

[1]

Дерфлингер, Герхард, Вольфганг Хёрманн и Йозеф Лейдольд. «Генерация случайных величин путем численной инверсии, когда известна только плотность». ACM Transactions on Modeling and Computer Simulation (TOMACS) 20.4 (2010): 1-25.