DensityDist from scipy.stats distribution

Thanks!
That is certainly an interesting possible way. I gave this a try and at least got the model to run but I must be doing something wrong as the fit is really poor and scipy was able to get a good fit.

Here is my pymc3 part:

def logp(value):
     return value - pm.math.exp(value)

with pm.Model() as model:
    scale = pm.HalfNormal('scale', sd=10)
    loc = pm.Uniform('loc', -100, 100)
    y = (dst.loc['2000'].values - loc) / scale
    gumbel_l = pm.DensityDist('gumbel_l', logp, testval=0, observed=y)
    trace = pm.sample(2000, )

Dealing with the loc and scale is a bit funny, is this done correctly? (See docs) I think I do have this wrong but getting errors other ways, see the notebook below.

The entire working example is at: