I ran in the meantime a MCMC with a Differential evolution algorithm (with sampling from the past but without snooker updates). I used initial values that were generated with `Metaheuristics.jl`

.

I got an acceptance rate between 5 and 10 % which is not too bad. The chains did not converge at least for some parameters (I ran two independent chains with three internal chains and a thinning rate of 10, red line on the right side is the prior). I have to check if reparametrisation makes sense:

I wanted to try the `Pigeons.jl`

package. But I got an error message:

```
ERROR: SliceSampler supports contrained target, but the sampler should be initialized in the support:
```

I think the problem is that not for all values that are possible from the priors a likelihood can be calculated and therefore I just return `-Inf`

in these cases. This happens for example if the overall biomass of a grassland is zero, then I cannot evaluate the likelihood of community weighted mean traits.

Is there a way to specify initial parameter values for the sampler?

Here is a reproducible example. I know it is not short and the precompilation takes quite some time.

```
import Pkg
Pkg.add(url="https://github.com/felixnoessler/GrasslandTraitSim.jl")
using Pigeons, Distributions, DistributionsAD, Turing
import GrasslandTraitSim.Valid as valid
import GrasslandTraitSim as sim
###################
training_plots = ["$(explo)$(lpad(i, 2, "0"))" for i in 1:9 for explo in ["HEG"]]
input_objs = valid.validation_input_plots(;
plotIDs = training_plots,
nspecies = 25,
trait_seed = 123);
valid_data = valid.get_validation_data_plots(;
plotIDs = training_plots);
preallocated = [sim.preallocate_vectors(; input_obj = input_objs[training_plots[1]])
for t in 1:Threads.nthreads()];
mp = valid.model_parameters(;)
###################
function turing_ll(x)
loglik_vec = Array{Float64}(undef, length(training_plots))
inf_p = (; zip(Symbol.(mp.names), x)...)
Threads.@threads for p in eachindex(training_plots)
loglik_vec[p] = valid.loglikelihood_model(sim; input_objs, valid_data,
calc = preallocated[Threads.threadid()],
plotID = training_plots[p], inf_p)
end
return sum(loglik_vec)
end
@model function calib_model()
moistureconv_alpha ~ Normal(0.0, 1.0)
moistureconv_beta ~ Normal(0.0, 1.0)
sen_α ~ truncated(Normal(0.0, 0.1); lower = 0)
sen_leaflifespan ~ truncated(Normal(0.0, 0.1); lower = 0)
sla_tr ~ truncated(Normal(0.02, 0.01); lower = 0)
sla_tr_exponent ~ truncated(Normal(1.0, 5.0); lower = 0)
βₚₑₜ ~ truncated(Normal(1.0, 1.0); lower = 0)
biomass_dens ~ truncated(Normal(1000.0, 1000.0); lower = 0)
belowground_density_effect ~ truncated(Normal(1.0, 0.5); lower = 0)
height_strength ~ Uniform(0.0, 1.0)
leafnitrogen_graz_exp ~ truncated(Normal(1.0, 5.0); lower = 0)
trampling_factor ~ truncated(Normal(200.0, 100.0); lower = 0)
grazing_half_factor ~ truncated(Normal(150.0, 50.0); lower = 0)
mowing_mid_days ~ truncated(Normal(10.0, 10.0); lower = 0)
totalN_β ~ truncated(Normal(0.1, 0.1); lower = 0)
CN_β ~ truncated(Normal(0.1, 0.1); lower = 0)
max_rsa_above_water_reduction ~ Uniform(0.0, 1.0)
max_sla_water_reduction ~ Uniform(0.0, 1.0)
max_amc_nut_reduction ~ Uniform(0.0, 1.0)
max_rsa_above_nut_reduction ~ Uniform(0.0, 1.0)
b_biomass ~ InverseGamma(5.0, 2000.0)
b_soilmoisture ~ truncated(Normal(0.0, 15.0); lower = 0)
b_sla ~ InverseGamma(10.0, 0.1)
b_lncm ~ InverseGamma(2.0, 5.0)
b_amc ~ InverseGamma(10.0, 3.0)
b_height ~ InverseGamma(10.0, 3.0)
b_rsa_above ~ InverseGamma(20.0, 0.2)
p = [moistureconv_alpha, moistureconv_beta, sen_α, sen_leaflifespan,
sla_tr, sla_tr_exponent, βₚₑₜ, biomass_dens, belowground_density_effect,
height_strength, leafnitrogen_graz_exp, trampling_factor,
grazing_half_factor, mowing_mid_days, totalN_β, CN_β,
max_rsa_above_water_reduction, max_sla_water_reduction,
max_amc_nut_reduction, max_rsa_above_nut_reduction,
b_biomass, b_soilmoisture, b_sla, b_lncm, b_amc, b_height, b_rsa_above]
Turing.@addlogprob! turing_ll(p)
return nothing
end
m = calib_model()
#### just check that the model can return a loglik
inf_p = (; zip(Symbol.(mp.names), mp.best)...)
TuringLogPotential(m)(inf_p)
#### use pigeons to run the model
pt = pigeons(target = TuringLogPotential(m))
```