Hi, I’m using latent GPs for a hierarchical extreme spatial model whose marginals are GEV distributions. So, I’m using also `pymc_experimental`

for GEV distributions.

My code is as follows.

```
with pm.Model() as model:
ell = pm.InverseGamma('ell', alpha=2, beta=1, shape=(3, 20))
# eta = pm.HalfNormal('eta', sigma=1, shape=3)
cov_shape = pm.gp.cov.Exponential(input_dim=20, ls=ell[0])
cov_loc = pm.gp.cov.Exponential(input_dim=20, ls=ell[1])
cov_scale = pm.gp.cov.Exponential(input_dim=20, ls=ell[2])
gp_shape = pm.gp.Latent(cov_func=cov_shape)
gp_loc = pm.gp.Latent(cov_func=cov_loc)
gp_log_scale = pm.gp.Latent(cov_func=cov_scale)
xi = gp_shape.prior('xi', X=X_)
mu = gp_loc.prior('mu', X=X_)
log_sigma = gp_log_scale.prior('log_sigma', X=X_)
# sigma = pm.Deterministic('sigma', pm.math.exp(log_sigma))
yy = pmx.distributions.GenExtreme('yy', mu=mu, sigma=pm.math.exp(log_sigma), xi=xi, observed=ys_maxima.T.ravel())
trace = pm.sample(1000, target_accept=0.95, random_seed=10)
```

Then, it returns the following error.

```
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Output exceeds the size limit. Open the full output data in a text editor---------------------------------------------------------------------------
SamplingError Traceback (most recent call last)
Cell In[74], line 2
1 with model:
----> 2 mp = pm.sample(1000, target_accept=0.95, random_seed=10)
File c:\XXX\Lib\site-packages\pymc\sampling\mcmc.py:619, in sample(draws, tune, chains, cores, random_seed, progressbar, step, nuts_sampler, initvals, init, jitter_max_retries, n_init, trace, discard_tuned_samples, compute_convergence_checks, keep_warning_stat, return_inferencedata, idata_kwargs, nuts_sampler_kwargs, callback, mp_ctx, model, **kwargs)
617 ip: Dict[str, np.ndarray]
618 for ip in initial_points:
--> 619 model.check_start_vals(ip)
620 _check_start_shape(model, ip)
622 # Create trace backends for each chain
File c:\XXX\Lib\site-packages\pymc\model.py:1779, in Model.check_start_vals(self, start)
1776 initial_eval = self.point_logps(point=elem)
1778 if not all(np.isfinite(v) for v in initial_eval.values()):
-> 1779 raise SamplingError(
1780 "Initial evaluation of model at starting point failed!\n"
1781 f"Starting values:\n{elem}\n\n"
1782 f"Logp initial evaluation results:\n{initial_eval}"
1783 )
SamplingError: Initial evaluation of model at starting point failed!
Starting values:
{'ell_log__': array([[ 0.01571654, -0.34856599, -0.49873131, -0.57287988, 0.95534123,
...
-0.89331143, 0.72452945, 0.37501074, -0.41410939, 0.51900205,
0.98783747, -0.82823462, -0.35654576, -0.91625861, -0.41124514])}
Logp initial evaluation results:
{'ell': -66.43, 'xi_rotated_': -1089.17, 'mu_rotated_': -1090.66, 'log_sigma_rotated_': -1084.57, 'yy': -inf}
```

Well, I took exponential to `log_sigma_rotated`

in different ways.

i) Defining `sigma = pm.Deterministic('sigma', pm.math.exp(log_sigma))`

;

ii) `sigma=pm.math.exp(log_sigma)`

and pass it to `GenExtreme`

’s sigma; and

iii) as shown in the above code.

But all of these return the same error.

FYI, I specified the random_seed in `pm.sample`

to try different initializations just in case. Interestingly, inference for the same hierarchical model for 1-dimensional simple simulation works.