CODEWITHSUNDEEP
bayesianpymc3
machadoyang
machadoyang

Reputation: 11

Sampling from posterior using custom likelihood in pymc3

I'm trying to create a custom likelihood using pymc3. The distribution is called Generalized Maximum Likelihood (GEV) which has the location (loc), scale (scale) and shape (c) parameters. The main ideia is to choose a beta distribution as a prior to the scale parameter and fix the location and scale parameters in the GEV likelihood. The GEV distribuition is not contained in the pymc3 standard distributions, so I have to create a custom likelihood. I googled it and found out that I should use the densitydist method but I don't know why it is incorrect.

See the code below:

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)

I'm sorry in advance if this is a dumb question, but I could'nt figure it out.

I'm studying annual maximum flows and I'm trying to implement the methodology described in "Generalized maximum-likelihood generalized extreme-value quantile estimators for hydrologic data" written by Martins and Stedinger.

Upvotes: 1

Views: 1464

Answers (1)

aseyboldt
aseyboldt

Reputation: 1090

If you mean the generalized extreme value distribution (https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution), then something like this should work (for 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)

Your logp function was invalid. Exponentiation is ** in python, and part of the expression wasn't valid for negative values.

Upvotes: 3

Related Questions