Вопрос по python, numpy, scipy – Создание новых дистрибутивов в scipy

11

Я пытаюсь создать дистрибутив на основе имеющихся у меня данных, а затем случайным образом извлечь из этого дистрибутива. Вот что у меня есть:

from scipy import stats
import numpy

def getDistribution(data):
    kernel = stats.gaussian_kde(data)
    class rv(stats.rv_continuous):
        def _cdf(self, x):
            return kernel.integrate_box_1d(-numpy.Inf, x)
    return rv()

if __name__ == "__main__":
    # pretend this is real data
    data = numpy.concatenate((numpy.random.normal(2,5,100), numpy.random.normal(25,5,100)))
    d = getDistribution(data)

    print d.rvs(size=100) # this usually fails

Я думаю, что это делает то, что я хочу, но я часто получаю ошибку (см. Ниже), когда я пытаюсь сделатьd.rvs(), а такжеd.rvs(100) никогда не работает. Я делаю что-то неправильно? Есть ли более простой или лучший способ сделать это? Если это ошибка в scipy, есть ли способ ее обойти?

Наконец, есть ли еще документация по созданию пользовательских дистрибутивов? Лучшее, что я нашел, - это документация scipy.stats.rv_continuous, которая довольно спартанская и не содержит полезных примеров.

След:

Traceback (most recent call last): File "testDistributions.py", line 19, in print d.rvs(size=100) File "/usr/local/lib/python2.6/dist-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/stats/distributions.py", line 696, in rvs vals = self._rvs(*args) File "/usr/local/lib/python2.6/dist-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/stats/distributions.py", line 1193, in _rvs Y = self._ppf(U,*args) File "/usr/local/lib/python2.6/dist-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/stats/distributions.py", line 1212, in _ppf return self.vecfunc(q,*args) File "/usr/local/lib/python2.6/dist-packages/numpy-1.6.1-py2.6-linux-x86_64.egg/numpy/lib/function_base.py", line 1862, in call theout = self.thefunc(*newargs) File "/usr/local/lib/python2.6/dist-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/stats/distributions.py", line 1158, in _ppf_single_call return optimize.brentq(self._ppf_to_solve, self.xa, self.xb, args=(q,)+args, xtol=self.xtol) File "/usr/local/lib/python2.6/dist-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/optimize/zeros.py", line 366, in brentq r = _zeros._brentq(f,a,b,xtol,maxiter,args,full_output,disp) ValueError: f(a) and f(b) must have different signs

Edit

Для тех, кому интересно, следуя совету в ответе ниже, вот код, который работает:

from scipy import stats
import numpy

def getDistribution(data):
    kernel = stats.gaussian_kde(data)
    class rv(stats.rv_continuous):
        def _rvs(self, *x, **y):
            # don't ask me why it's using self._size 
            # nor why I have to cast to int
            return kernel.resample(int(self._size)) 
        def _cdf(self, x):
            return kernel.integrate_box_1d(-numpy.Inf, x)
        def _pdf(self, x):
            return kernel.evaluate(x)
    return rv(name='kdedist', xa=-200, xb=200)
Поэтому, когда мы делаем вышеупомянутое иrandoms = getDistribution(Mydata) а потомrandoms = randoms.rvs(size=1000) это выполняет триdef внутри класса? Т.е. расчет pdf, его интеграция и т. д.? ThePredator
Я заставляю своих случайных людей следить за распределением данных, но я бы хотел сгладить их так, чтобы они не точно соответствовали распределению данных. Я вручную настраивал пропускную способность вkernel сделать это. Например, что-то вроде того, как мы указываем функцию PDF, а затем используем функцию PDF для создания случайных событий с помощью Metropolis Hastings. ThePredator

Ваш Ответ

1   ответ
7

rvs использует инверсию cdf, ppf, для создания случайных чисел. Так как вы не указываете ppf, он рассчитывается по алгоритму rootfinding,brentq. brentq использует нижнюю и верхнюю границы того, где следует искать значение, при котором функция равна нулю (найдите x такой, что cdf (x) = q, q - квантиль).

По умолчанию для лимитов,xa а такжеxbслишком малы в вашем примере. Следующее работает для меня со Scipy 0.9.0,xa, xb может быть установлен при создании экземпляра функции

def getDistribution(data):
    kernel = stats.gaussian_kde(data)
    class rv(stats.rv_continuous):
        def _cdf(self, x):
            return kernel.integrate_box_1d(-numpy.Inf, x)
    return rv(name='kdedist', xa=-200, xb=200)

В настоящее время есть запрос на извлечение scipy, чтобы улучшить это, поэтому в следующем выпускеxa а такжеxb будет расширен автоматически, чтобы избежатьf(a) and f(b) must have different signs исключение.

По этому вопросу не так много документации, проще всего следовать некоторым примерам (и спрашивать в списке рассылки).

редактировать: дополнение

pdf: Поскольку у вас есть функция плотности, также заданная gaussian_kde, я бы добавил_pdf метод, который сделает некоторые расчеты более эффективными.

edit2: дополнение

rvs: Если вы заинтересованы в генерации случайных чисел, то у gaussian_kde есть метод повторной выборки. Случайная выборка может быть сгенерирована путем выборки из данных и добавления гауссовского шума. Таким образом, это будет быстрее, чем общие rvs, использующие метод ppf. Я написал бы метод ._rvs, который просто вызывает метод resample gaussian_kde.

precomputing ppfЯ не знаю ни одного общего способа предварительного вычисления PPF. Однако способ, которым я думал об этом (но никогда не пытался), состоит в том, чтобы предварительно вычислить ppf во многих точках, а затем использовать линейную интерполяцию для приближения функции ppf.

edit3: о_rvs ответить на вопрос Шриватсана в комментарии

_rvs является специфичным для распределения методом, который вызывается открытым методомrvs. rvs является универсальным методом, который выполняет некоторую проверку аргументов, добавляет местоположение и масштаб, и устанавливает атрибутself._size который является размером запрошенного массива случайных переменных, а затем вызывает метод, специфичный для распределения._rvs или это общий аналог. Дополнительные аргументы в._rvs являются параметрами формы, но так как в этом случае их нет,*x а также**y избыточны и не используются.

Я не знаю, насколько хорошоsize или форма.rvs Метод работает в многомерном случае. Эти дистрибутивы предназначены для одномерных дистрибутивов и могут не полностью работать для многовариантного случая, или могут потребоваться некоторые изменения.

Error: User Rate Limit Exceeded
Error: User Rate Limit Exceeded Noah
Error: User Rate Limit Exceeded_rvs
Error: User Rate Limit Exceeded
Error: User Rate Limit Exceeded Noah

Похожие вопросы