if (FALSE) {
  # Run this part of the code only once
  install.packages("embryogrowth")
  
  install.packages("https://hebergement.universite-paris-saclay.fr/marcgirondot/CRAN/HelpersMG.tar.gz", repos=NULL, type="source")
  install.packages("https://hebergement.universite-paris-saclay.fr/marcgirondot/CRAN/embryogrowth.tar.gz", repos=NULL, type="source")
}

library(embryogrowth)
Lepidochelys.olivacea = subset(DatabaseTSD, Species=="Lepidochelys olivacea" & (!is.na(Sexed) & Sexed!=0))
colnames(Lepidochelys.olivacea)
unique(Lepidochelys.olivacea$RMU.2023)
nrow(Lepidochelys.olivacea)

Lepidochelys.olivacea.WA = subset(Lepidochelys.olivacea, subset=RMU.2023 == "West Atlantic")
nrow(Lepidochelys.olivacea.WA)

# Here the intersexes are excluded from the analysis

tsdF = tsd(df=Lepidochelys.olivacea.WA, 
           equation = "flexit**", 
           parameters.initial = c(P=30, SH=1, SL=1))
plot(tsdF)

priorsN = tsd_MHmcmc_p(result = tsdF, default = "dnorm", accept = TRUE)
priorsU = tsd_MHmcmc_p(result = tsdF, default = "dunif", accept = TRUE)
priors = rbind(priorsN[1, , drop=FALSE], priorsU[2:3, ])
priors[1, "Prior1"] = 30
priors[2:3, c("Prior1", "Min")] = 0
priors[2:3, c("Prior2", "Max")] = 10
priors

tsdF_mcmc = tsd_MHmcmc(result = tsdF, parametersMCMC = priors, 
           n.iter = 20000, thin = 10, n.adapt=1000)

plot(tsdF_mcmc, parameters = "P", what="MarkovChain")
plot(tsdF_mcmc, parameters = "SL", what="MarkovChain", main="SL")
plot(tsdF_mcmc, parameters = "SH", what="MarkovChain", main="SH")

as.parameters(x=tsdF_mcmc, index = "quantile", probs = c(0.025, 0.5, 0.975))
as.parameters(x=tsdF_mcmc, index = "quantile", probs = c(0.25, 0.5, 0.75))
plot(tsdF, resultmcmc = tsdF_mcmc, replicate.CI = 2000)