if (FALSE) {
library(phenology)
# Example
data <- c(0, 47, 15, 6, 5, 4, 2, 5, 57, 203, 205, 103, 35, 24, 12, 10,
13, 49, 86, 107, 111, 73, 47, 30, 19, 17, 33, 48, 77, 83, 65,
37, 27, 23, 24, 22, 41, 42, 44, 33, 39, 24, 18, 18, 22, 22, 19,
24, 28, 17, 18, 19, 17, 4, 12, 9, 6, 11, 7, 11, 12, 5, 4, 6,
11, 5, 6, 7, 3, 2, 1, 3, 2, 1, 2, 0, 0, 3, 1, 0, 2, 0, 0, 1)
class(data) <- unique(append("IP", class(data)))
plot(data)
######### Fit parametric ECF using Maximum-Likelihood
par <- c(meanIP = 9.9959691992722917,
sdIP = 0.10066664270893474,
minIP = 7.5684588178888754,
pAbort = 2.2510012544630911,
meanAbort = 2.8969185085603386,
sdAbort = 0.92688983853803242,
pCapture = -1.0393803705929086,
meanECF = 3.9551519427394255,
sdECF = 0.31657679943365019)
fML <- IPFit(x=par,
fixed.parameters=c(N=1000000),
data=data,
verbose=FALSE,
model="ML")
# Plot the fitted ECF
plot(fML, result="ECF")
# Plot the Internesting Period distribution
plot(fML, result="IP")
# Plot the distribution of days between tentatives
plot(fML, result="Abort")
plot(fML, result="Abort", xlim=c(0, 10))
# Plot the data
plot(fML, result="data")
# Plot the data and the model
plot(fML, result="data&model")
# Plot the cumulative proportion of ECF according to date of observation
plot(fML, result="reverseECF")
plot(fML_NP_Delta, result="reverseECF", col=grey.colors(128))
######### Fit using Metropolis-Hastings
# ECF.1 = 1 is fixed
par <- c(ECF.2 = 0.044151921569961131,
ECF.3 = 2.0020778325280748,
ECF.4 = 2.6128345101617083,
ECF.5 = 2.6450582416622375,
ECF.6 = 2.715198206774927,
ECF.7 = 2.0288031327239904,
ECF.8 = 1.0028041546528881,
ECF.9 = 0.70977432157689235,
ECF.10 = 0.086052204035003091,
ECF.11 = 0.011400419961702518,
ECF.12 = 0.001825219438794076,
ECF.13 = 0.00029398731859899116,
ECF.14 = 0.002784886479846703,
meanIP = 9.9887100433529721,
sdIP = 0.10580250625108811,
minIP = 6.5159124624132048,
pAbort = 2.5702251748938956,
meanAbort = 2.2721679285648841,
sdAbort = 0.52006431730489933,
pCapture = 0.079471782729506113)
df <- data.frame(Density=rep("dunif", length(par)),
Prior1=c(rep(0, 13), 8, 0.001, 0, -8, 0, 0.001, -8),
Prior2=c(rep(10, 13), 12, 1, 10, 8, 2, 1, 8),
SDProp=unname(c(ECF.2 = 6.366805760909012e-05,
ECF.3 = 6.366805760909012e-05,
ECF.4 = 6.366805760909012e-05,
ECF.5 = 6.366805760909012e-05,
ECF.6 = 6.366805760909012e-05,
ECF.7 = 6.366805760909012e-05,
ECF.8 = 6.366805760909012e-05,
ECF.9 = 6.366805760909012e-05,
ECF.10 = 6.366805760909012e-05,
ECF.11 = 6.366805760909012e-05,
ECF.12 = 6.366805760909012e-05,
ECF.13 = 6.366805760909012e-05,
ECF.14 = 6.366805760909012e-05,
meanIP = 6.366805760909012e-05,
sdIP = 6.366805760909012e-05,
minIP = 6.366805760909012e-05,
pAbort = 6.366805760909012e-05,
meanAbort = 6.366805760909012e-05,
sdAbort = 6.366805760909012e-05,
pCapture = 6.366805760909012e-05)),
Min=c(rep(0, 13), 8, 0.001, 0, -8, 0, 0.001, -8),
Max=c(rep(10, 13), 12, 1, 10, 8, 2, 1, 8),
Init=par, stringsAsFactors = FALSE)
rownames(df)<- names(par)
fMH <- IPFit(parametersMH=df,
fixed.parameters=c(N=10000),
data=data,
verbose=FALSE,
n.iter = 10000,
n.chains = 1, n.adapt = 100, thin = 1, trace = TRUE,
adaptive = TRUE,
model="MH")
# Plot the fitted ECF
plot(fMH, result="ECF")
# Plot the posteriors and priors
plot(fMH$MH, parameters="meanIP", xlim=c(6, 14))
plot(x=1:length(fMH$MH$resultLnL[[1]]), y=fMH$MH$resultLnL[[1]],
type="l", xlab="Iterations", ylab="Ln L", bty="n", las=1)
}
Run the code above in your browser using DataLab