Выборка из задней части с использованием пользовательского правдоподобия в pymc3

Question

Я пытаюсь создать пользовательскую вероятность, используя pymc3. Распределение называется Обобщенным максимальным правдоподобием (GEV), которое имеет параметры местоположения (loc), масштаба (scale) и shape (c). Основная идея состоит в том, чтобы выбрать бета-распределение в качестве предшествующего параметра масштаба и зафиксировать параметры местоположения и масштаба в вероятности GEV. Дистрибутив GEV не содержится в стандартных дистрибутивах pymc3, поэтому я должен создать собственную вероятность. Я прогуглил это и обнаружил, что должен использовать метод densitydist, но я не знаю, почему это неверно.

Смотрите код ниже:

import pymc3 as pm
import numpy as np
from theano.tensor import exp

data=np.random.randn(20)

with pm.Model() as model:
    c=pm.Beta('c',alpha=6,beta=9)
    loc=1
    scale=2
    gev=pm.DensityDist('gev', lambda value: exp(-1+c*(((value-loc)/scale)^(1/c))), testval=1)
    modelo=pm.gev(loc=loc, scale=scale, c=c, observed=data)
    step = pm.Metropolis()
    trace = pm.sample(1000, step)
pm.traceplot(trace)

Я заранее извиняюсь, если это глупый вопрос, но я не могу понять это.

Я изучаю годовые максимальные потоки и пытаюсь реализовать методологию, описанную в "Обобщенных квантильных оценках экстремальных значений максимального правдоподобия для гидрологических данных", написанных Мартинсом и Стедингером.

Answer 1

Если вы имеете в виду обобщенное распределение экстремальных значений ( https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution), то что-то вроде этого должно работать (для c!= 0):

import pymc3 as pm
import numpy as np
import theano.tensor as tt
from pymc3.distributions.dist_math import bound

data = np.random.randn(20)


with pm.Model() as model:
    c = pm.Beta('c', alpha=6, beta=9)
    loc = 1
    scale = 2

    def gev_logp(value):
        scaled = (value - loc) / scale
        logp = -(scale
                 + ((c + 1) / c) * tt.log1p(c * scaled)
                 + (1 + c * scaled) ** (-1/c))
        alpha = loc - scale / c
        bounds = tt.switch(value > 0, value > alpha, value < alpha)
        return bound(logp, bounds, c != 0)

    gev = pm.DensityDist('gev', gev_logp, observed=data)
    trace = pm.sample(2000, tune=1000, njobs=4)
pm.traceplot(trace)

Ваша функция logp была недействительной. Экспонирование ** в Python, а часть выражения не действительны для отрицательных значений.