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))

###########################################
# Here the intersexes are included as males
###########################################
males <- Lepidochelys.olivacea$Males + ifelse(is.na(Lepidochelys.olivacea$Intersexes), 0, Lepidochelys.olivacea$Intersexes)
females <- Lepidochelys.olivacea$Females
temperatures <- Lepidochelys.olivacea$Incubation.temperature.set

tsdF_Intersexes <- tsd(males = males                                                       , 
            females = females                                                   , 
            temperatures = temperatures                                         , 
            equation = "flexit**", 
            parameters.initial = c(P=30, SH=1, SL=1))

priorsN = tsd_MHmcmc_p(result = tsdF_Intersexes, default = "dnorm", accept = TRUE)
priorsU = tsd_MHmcmc_p(result = tsdF_Intersexes, 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

tsdF_Intersexes_mcmc = tsd_MHmcmc(result = tsdF_Intersexes, parametersMCMC = priors, 
                       n.iter = 20000, thin = 10, n.adapt=1000)
plot(tsdF_Intersexes, resultmcmc = tsdF_Intersexes_mcmc, replicate.CI = 2000, xlim=c(24, 36))

####################################################
# Here the intersexes are excluded from the analysis
####################################################

males <- Lepidochelys.olivacea$Males
females <- Lepidochelys.olivacea$Females
temperatures <- Lepidochelys.olivacea$Incubation.temperature.set

tsdF_NoIntersexes <- tsd(males = males                                                       , 
            females = females                                                   , 
            temperatures = temperatures                                         , 
            equation = "flexit**", 
            parameters.initial = c(P=30, SH=1, SL=1))

priorsN = tsd_MHmcmc_p(result = tsdF_NoIntersexes, default = "dnorm", accept = TRUE)
priorsU = tsd_MHmcmc_p(result = tsdF_NoIntersexes, 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

tsdF_NoIntersexes_mcmc = tsd_MHmcmc(result = tsdF_NoIntersexes, parametersMCMC = priors, 
                                  n.iter = 20000, thin = 10, n.adapt=1000)
plot(tsdF_NoIntersexes, resultmcmc = tsdF_NoIntersexes_mcmc, replicate.CI = 2000, xlim=c(24, 36))