Преобразованное отклонение плотности (TDR)

• Требуется: T-вогнутая PDF, dPDF

• Опционально: mode, center

• Скорость:

  • Настройка: медленно

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

TDR — это метод принятия/отклонения, который использует вогнутость преобразованной плотности для построения функции-шапки и автоматического сжатия. Такие PDF называются T-вогнутыми. В настоящее время реализованы следующие преобразования:

\[\begin{split}c = 0.: T(x) &= \log(x)\\ c = -0.5: T(x) &= \frac{1}{\sqrt{x}} \text{ (Default)}\end{split}\]

В дополнение к PDF, также требуется производная PDF по x (т.е. вариант). Эти функции должны присутствовать как методы объекта Python, который затем можно передать генераторам для создания их объекта. Реализованный вариант использует сжатия, пропорциональные функции-шляпе ([1]).

Пример использования этого метода показан ниже:

>>> import numpy as np
>>> from scipy.stats.sampling import TransformedDensityRejection
>>> from scipy.stats import norm
>>>
>>> class StandardNormal:
...     def pdf(self, x):
...         # note that the normalization constant is not required
...         return np.exp(-0.5 * x*x)
...     def dpdf(self, x):
...         return -x * np.exp(-0.5 * x*x)
...
>>> dist = StandardNormal()
>>>
>>> urng = np.random.default_rng()
>>> rng = TransformedDensityRejection(dist, random_state=urng)
>>> rng.rvs()
-1.526829048388144

В приведённом выше примере мы использовали метод TDR для выборки из стандартного нормального распределения. Обратите внимание, что мы можем опустить константу нормализации при вычислении PDF. Обычно это помогает ускорить этап выборки. Также обратите внимание, что PDF не нужно векторизовать. Он должен принимать и возвращать скаляр.

Также возможно напрямую вычислить обратную функцию CDF распределения hat используя ppf_hat метод.

>>> rng.ppf_hat(0.5)
-0.00018050266342362759
>>> norm.ppf(0.5)
0.0
>>> u = np.linspace(0, 1, num=10)
>>> rng.ppf_hat(u)
array([       -inf, -1.22227372, -0.7656556 , -0.43135731, -0.14002921,
        0.13966423,  0.43096141,  0.76517113,  1.22185606,         inf])
>>> norm.ppf(u)
array([       -inf, -1.22064035, -0.76470967, -0.4307273 , -0.1397103 ,
        0.1397103 ,  0.4307273 ,  0.76470967,  1.22064035,         inf])

Помимо метода PPF, другие атрибуты могут быть доступны чтобы увидеть, насколько хорошо генератор соответствует заданному распределению. Это:

• ‘squeeze_hat_ratio’: (площадь под squeeze) / (площадь под hat) для генератора. Это число между 0 и 1. Ближе к 1 означает, что функции hat и squeeze плотно охватывают распределение, и требуется меньше вычислений PDF для генерации выборок. Ожидаемое количество вычислений плотности ограничено (1/squeeze_hat_ratio) - 1 на выборку. По умолчанию он поддерживается выше 0.99, но это можно изменить, передав max_squeeze_hat_ratio параметр.

• ‘hat_area’: площадь под шляпой для генератора.

• ‘squeeze_area’: область под сжатием для генератора.

  >>> rng.squeeze_hat_ratio
  0.9947024204884917
  >>> rng.hat_area
  2.510253139791547
  >>> rng.squeeze_area
  2.4969548741894876
  >>> rng.squeeze_hat_ratio == rng.squeeze_area / rng.hat_area
  True
  

Распределение может быть усечено путем передачи параметра domain:

>>> urng = np.random.default_rng()
>>> rng = TransformedDensityRejection(dist, domain=[0, 1], random_state=urng)
>>> rng.rvs(10)
array([0.05452512, 0.97251362, 0.49955877, 0.82789729, 0.33048885,
       0.55558548, 0.23168323, 0.13423275, 0.73176575, 0.35739799])

Если домен не указан, support метод dist объект используется для определения области:

>>> class StandardNormal:
...     def pdf(self, x):
...         return np.exp(-0.5 * x*x)
...     def dpdf(self, x):
...         return -x * np.exp(-0.5 * x*x)
...     def support(self):
...         return -np.inf, np.inf
...
>>> dist = StandardNormal()
>>>
>>> urng = np.random.default_rng()
>>> rng = TransformedDensityRejection(dist, random_state=urng)
>>> rng.rvs(10)
array([-1.52682905,  2.06206883,  0.15205036,  1.11587367, -0.30775562,
       0.29879802, -0.61858268, -1.01049115,  0.78853694, -0.23060766])

Если dist объект не предоставляет support метод, предполагается, что область определения (-np.inf, np.inf).

Для увеличения squeeze_hat_ratio, передать max_squeeze_hat_ratio:

>>> dist = StandardNormal()
>>> rng = TransformedDensityRejection(dist, max_squeeze_hat_ratio=0.999,
...                                   random_state=urng)
>>> rng.squeeze_hat_ratio
0.999364900465214

Давайте посмотрим, как это влияет на обратные вызовы к методу PDF распределения:

>>> from copy import copy
>>> class StandardNormal:
...     def __init__(self):
...         self.callbacks = 0
...     def pdf(self, x):
...         self.callbacks += 1
...         return np.exp(-0.5 * x*x)
...     def dpdf(self, x):
...         return -x * np.exp(-0.5 * x*x)
...
>>> dist1 = StandardNormal()
>>> urng1 = np.random.default_rng()
>>> urng2 = copy(urng1)
>>> rng1 = TransformedDensityRejection(dist1, random_state=urng1)
>>> dist1.callbacks  # evaluations during setup
139
>>> dist1.callbacks = 0  # don't consider evaluations during setup
>>> rvs = rng1.rvs(100000)
>>> dist1.callbacks  # evaluations during sampling
527
>>> dist2 = StandardNormal()
>>> # use the same stream of uniform random numbers
>>> rng2 = TransformedDensityRejection(dist2, max_squeeze_hat_ratio=0.999,
...                                    random_state=urng2)
>>> dist2.callbacks  # evaluations during setup
467
>>> dist2.callbacks = 0  # don't consider evaluations during setup
>>> rvs = rng2.rvs(100000)
>>> dist2.callbacks  # evaluations during sampling
84      # may vary

Как мы видим, при увеличении требуется значительно меньше вычислений PDF во время выборки. squeeze_hat_ratio. Функция PPF-hat также более точна:

>>> abs(norm.ppf(0.975) - rng1.ppf_hat(0.975))
0.0027054565421578136
>>> abs(norm.ppf(0.975) - rng2.ppf_hat(0.975))
0.00047824084476300044

Однако обратите внимание, что это происходит за счёт увеличения вычислений PDF во время настройки.

Для плотностей с модами, не близкими к 0, рекомендуется задать либо моду, либо центр распределения, передав mode или center параметры. Последнее — приблизительное местоположение моды или среднего распределения. Это местоположение предоставляет некоторую информацию об основной части PDF и используется для избежания численных проблем.

>>> # mode = 0 for our distribution
>>> # if exact mode is not available, pass 'center' parameter instead
>>> rng = TransformedDensityRejection(dist, mode=0.)

По умолчанию метод использует 30 точек построения для создания шапки. Это можно изменить, передав construction_points параметр, который может быть либо массивом точек построения, либо целым числом, представляющим количество используемых точек построения.

>>> rng = TransformedDensityRejection(dist,
...                                   construction_points=[-5, 0, 5])

Этот метод принимает множество других параметров настройки. См. документацию для полного списка. Дополнительная информация о параметрах и методе может быть найдена в Раздел 5.3.16 руководства пользователя UNU.RAN.

Пожалуйста, смотрите [1] и [2] для получения дополнительных сведений об этом методе.

Ссылки

[1] (1,2)

Руководство по UNU.RAN, Раздел 5.3.16, 'TDR - Transformed Density Rejection', http://statmath.wu.ac.at/software/unuran/doc/unuran.html#TDR

[2]

Хёрманн, Вольфганг. «Техника отклонения для выборки из T-вогнутых распределений.» ACM Transactions on Mathematical Software (TOMS) 21.2 (1995): 182-193