I am trying to train a generalized extreme value model with the code below. The model is described in Stuart Coles book “An Introduction to Statistical Modeling of Extreme Values” and I have used the same data port pirie sea levels. The code is:

```
with pm.Model() as model_gev:
loc = pm.Normal('loc', mu=3, sigma=10)
sig=pm.Normal('sigma', mu=0.24, sigma=10)
xi = pm.Normal('xi', mu=0, sigma=10)
def gev_logp(value):
scaled = (value - loc) / sig
logp = -(log(sig)
+ ((xi + 1) / xi) * tt.log1p(xi * scaled)
+ (1 + xi * scaled) ** (-1/xi))
alpha = loc - sig/xi
bounds = tt.switch(xi > 0, xi > alpha, xi < alpha)
return bound(logp, bounds, xi!= 0)
gev = pm.DensityDist('gev', gev_logp, observed=series)
step = pm.Metropolis()
trace = pm.sample(10000, chains=4, tune=3000,step=step,start={'loc': 1.0, 'sigma': 1.0, 'xi': 1.0})
pm.traceplot(trace)
```

I get the following estimates which is very close to maximum likelihood estimates:

loc=3.88, sigma=0.21 and xi=-0.07. rhat is 1.

My question is the density plots of the parameters. Do they look reasonable?