Существует ли быстрая альтернатива scipy _norm_pdf для выборки с коррелированным распределением? - PullRequest
Новая
       38

Существует ли быстрая альтернатива scipy _norm_pdf для выборки с коррелированным распределением?

5 голосов
/ 01 февраля 2020

Я выбрал серию непрерывных распределений SciPy для моделирования Монте-Карло и собираюсь взять большое количество выборок из этих распределений. Однако я хотел бы иметь возможность брать коррелированные выборки, так что i -ая выборка берет, например, 90-й процентиль от каждого из распределений.

При этом я обнаружил причуду в производительности SciPy: 

# very fast way to many uncorrelated samples of length n
for shape, loc, scale, in distro_props:
    sp.stats.norm.rvs(*shape, loc=loc, scale=scale, size=n)

# verrrrryyyyy slow way to take correlated samples of length n
correlate = np.random.uniform(size=n)
for shape, loc, scale, in distro_props:
    sp.stats.norm.ppf(correlate, *shape, loc=loc, scale=scale)

Большинство результатов об этом утверждают, что медлительность в этих дистрибутивах SciPy из проверки типа и т. Д. c. оберток. Однако когда я профилировал код, большая часть времени была потрачена на основную математическую функцию [_continuous_distns.py:179(_norm_pdf)] 1 . Кроме того, он масштабируется с n, подразумевая, что он циклически проходит через каждый элемент.

Документы SciPy на rv_continuous почти предполагают, что подкласс должен переопределить это для производительности, но это кажется странным, что я бы запрятал в SciPy, чтобы ускорить их PPF. Я бы просто вычислил это для нормалей по формуле ppf, но я также использую lognormal и skewed normal, которые более трудны для реализации.

Итак, что является лучшим способом в Python для вычисления быстрый ppf для нормального, логнормального и искаженного нормального распределения? Или, в более широком смысле, взять коррелированные выборки из нескольких таких распределений?

Ответы [ 2 ]

3 голосов
/ 06 февраля 2020

Если вам нужен только обычный ppf, действительно удивительно, что он такой медленный, но вместо этого вы можете использовать scipy.special.erfinv: 

x = np.random.uniform(0,1,100)
np.allclose(special.erfinv(2*x-1)*np.sqrt(2),stats.norm().ppf(x))
# True
timeit(lambda:stats.norm().ppf(x),number=1000)
# 0.7717257660115138
timeit(lambda:special.erfinv(2*x-1)*np.sqrt(2),number=1000)
# 0.015020604943856597

РЕДАКТИРОВАТЬ:

lognormal и triangle также прямолинейны: 

c = np.random.uniform()

np.allclose(np.exp(c*special.erfinv(2*x-1)*np.sqrt(2)),stats.lognorm(c).ppf(x))
# True

np.allclose(((1-np.sqrt(1-(x-c)/((x>c)-c)))*((x>c)-c))+c,stats.triang(c).ppf(x))
# True

нормальный перекос Я, к сожалению, недостаточно знаком.

0 голосов
/ 03 марта 2020

В конечном итоге эта проблема была вызвана моим использованием дистрибутива skew-normal . У ppf косо-нормального на самом деле нет аналитического определения c в замкнутой форме, поэтому для вычисления ppf он прибегнул к приближению scipy.continuous_rv цифри c, что потребовало многократного вычисления cdf и используя это, чтобы сосредоточиться на значении PPF. Отклонение от нормального pdf является произведением нормального pdf и нормального cdf, поэтому это приближение чисел c многократно называлось pdf и cdf нормали. Вот почему, когда я профилировал код, он выглядел как проблема с нормальным распределением, а не с SKU. Другой ответ на этот вопрос смог добиться экономии времени за счет пропуска проверки типов, но на самом деле не имел значения для роста времени выполнения, а только для времени выполнения с малым n.

Чтобы решить эту проблему, я заменил асимметричный дистрибутив на дистрибутив Johnson SU . У него на 2 свободных параметра больше, чем у нормального распределения, поэтому он может эффективно соответствовать различным типам перекоса и эксцесса. Он поддерживается для всех действительных чисел и имеет определение ppf в закрытой форме с быстрой реализацией в SciPy. Ниже вы можете увидеть примеры дистрибутивов Johnson SU, которые я подбирал с 10-го, 50-го и 90-го процентилей.

Example distributions

...