Basic masting model
library(MASS)
library(rjags)Loading required package: coda
Linked to JAGS 4.3.0
Loaded modules: basemod,bugsSimulation
Relationship between the number of fruits at year n-1 and the probability of making buds at year n
We use a logistic function and reflect it to make it decrease with increasing fruits at year n-1.
# Logistic function: y goes from 0 to maximum as x increases
logistic <- function(x, midpoint, maximum=1, steepness=0.1) {
y = maximum/(1+exp(-steepness*(x-midpoint)))
return(y)
}
# Reflected logistic function: y goes from maximum to 0 as x increases
ref_logistic <- function(x, midpoint, maximum=1, steepness=0.1) {
return(maximum - logistic(x, midpoint, maximum, steepness))
}
Plotting the logistic function we use
x = (1:1000)*2
plot(x, ref_logistic(x, midpoint=500), t="l") 
# Simulates fructification at a given year with a simple model that depends on the previous year's number of fruits only.
# This dependency is modelled through the ref_logistic function, which provides a probability of making buds p_buds.
simulateFructification <- function (fruits_last_year, midpoint, maximum=1, steepness=0.1, meanFruits=1000, theta_fruits=10) {
p_buds = ref_logistic(fruits_last_year, midpoint, maximum, steepness)
#print(paste("Probability of getting buds this year: ", p_buds, sep=""))
n_fruits = rnegbin(1, mu=meanFruits * p_buds, theta=theta_fruits)
return(n_fruits)
}
# Simulates fructifications over N years, by calling iteratively simulateFructification.
simulateNYearsFructification <- function (Nyears=7, midpoint=500, maximum=1, steepness=0.1, meanFruits=1000, theta_fruits=10) {
fruits = c()
# we arbitrarily start assuming there had been fruits the year preceding the samples
fruits = c(fruits, rnegbin(1, mu=meanFruits, theta=theta_fruits))
for (year in 1:Nyears) {
fruits = c(fruits, simulateFructification(fruits[length(fruits)], midpoint, maximum, steepness, meanFruits, theta_fruits))
}
return(fruits)
}
# Simulates fructifications over N years and in several areas, by calling iteratively simulateNYearsFructification.
simulateNAreas <- function(Nareas=15) {
all_fruits = c()
for (area in 1:Nareas) {
all_fruits = rbind(all_fruits, simulateNYearsFructification())
}
return(all_fruits)
}
Simulating one area
fruits = simulateNYearsFructification()What does it look like?
plot(fruits, xlab="Year")
Simulating 15 areas
fruits = simulateNAreas()What does it look like?
plot(x=1:8, y=fruits[1,], xlab="Year", ylab="Number of fruits", t="l", col=0, ylim=c(0,2000))
for (i in 2:dim(fruits)[1]) {
lines(fruits[i,], col=i)
}
Estimation
We are going to use Jags to estimate the parameters of a model that matches the simulation model. The distribution dnegbin in Jags is parametrized in p and s, whereas R’s rnegbin is parametrized in mu (the mean) and theta. The relationships between those parameters are as follows: - s = theta - mu = s(1-p)/p i.e. p = s/(mu+s) In the simulation, we have: mu = p_buds * p1 * meanFruits We want to have these parameters in the Jags estimation model to be able to compare the estimated to the true values: the estimation model has to match the simulation model.
jags_model <-
"model
{
for (area in 1:Nareas)
{
# special case: year 1
fruits[area, 1] ~ dnegbin(min(scalePrior/(mu[area, 1]+scalePrior), 0.999), scalePrior)
mu[area, 1] <- pBuds[area, 1] * meanFruits
pBuds[area, 1] ~ dbeta(1.0,1.0)
# years 2 and following
for (year in 2:Nyears)
{
fruits[area, year] ~ dnegbin(min(scalePrior/(mu[area, year]+scalePrior), 0.999), scalePrior)
mu[area,year] <- pBuds[area, year] * meanFruits
pBuds[area, year] <- 1.0 - 1.0/(1+exp(-steepness*(fruits[area, year-1]-midpoint)))
}
}
scalePrior ~ dunif(0.1, 5000)
meanFruits ~ dunif(5, 5000)
midpoint ~ dunif(0, 5000)
steepness ~ dunif(0, 5)
}"dcomplete = list(fruits=fruits, Nareas=dim(fruits)[1], Nyears=dim(fruits)[2])
ini <- list(list(scalePrior=1, meanFruits = 500, midpoint = 500, steepness = 0.1),
list(scalePrior=10, meanFruits = 1000, midpoint = 1000, steepness = 1.0),
list(scalePrior=100, meanFruits = 100, midpoint = 100, steepness = 4.0))
m1<-jags.model(file = textConnection(jags_model), data = dcomplete, inits = ini, n.chains = 3)Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 120
Unobserved stochastic nodes: 19
Total graph size: 808
Initializing model
|
| | 0%
|
|+ | 2%
|
|++ | 4%
|
|+++ | 6%
|
|++++ | 8%
|
|+++++ | 10%
|
|++++++ | 12%
|
|+++++++ | 14%
|
|++++++++ | 16%
|
|+++++++++ | 18%
|
|++++++++++ | 20%
|
|+++++++++++ | 22%
|
|++++++++++++ | 24%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 28%
|
|+++++++++++++++ | 30%
|
|++++++++++++++++ | 32%
|
|+++++++++++++++++ | 34%
|
|++++++++++++++++++ | 36%
|
|+++++++++++++++++++ | 38%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 42%
|
|++++++++++++++++++++++ | 44%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 48%
|
|+++++++++++++++++++++++++ | 50%
|
|++++++++++++++++++++++++++ | 52%
|
|+++++++++++++++++++++++++++ | 54%
|
|++++++++++++++++++++++++++++ | 56%
|
|+++++++++++++++++++++++++++++ | 58%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 62%
|
|++++++++++++++++++++++++++++++++ | 64%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 68%
|
|+++++++++++++++++++++++++++++++++++ | 70%
|
|++++++++++++++++++++++++++++++++++++ | 72%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|++++++++++++++++++++++++++++++++++++++ | 76%
|
|+++++++++++++++++++++++++++++++++++++++ | 78%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 82%
|
|++++++++++++++++++++++++++++++++++++++++++ | 84%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 88%
|
|+++++++++++++++++++++++++++++++++++++++++++++ | 90%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 92%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|++++++++++++++++++++++++++++++++++++++++++++++++ | 96%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 98%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
# burnin
update(m1, 5000)
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
# sampling
mcmc1 <- coda.samples(m1, c("meanFruits", "midpoint", "steepness", "scalePrior"), n.iter = 10000, thin = 5)
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%Convergence diagnostics
gelman.diag(mcmc1)Potential scale reduction factors:
Point est. Upper C.I.
meanFruits 1 1
midpoint 1 1
scalePrior 1 1
steepness 1 1
Multivariate psrf
1The potential scale reduction factors seem close to 1.0 for both the mean and the upper C.I., which is a good sign that all 3 chains converged to similar posterior densities.
effectiveSize(mcmc1)meanFruits midpoint scalePrior steepness
5868.638 2271.160 5825.108 2137.073 It would seem like we have enough independent samples to estimate the posterior distributions well.
Analysis of the estimated parameters
plot(mcmc1)
We don’t see much, we’ll look at each variable in turn.
meanFruits
plot(mcmc1[,c("meanFruits")])
summary(mcmc1[,c("meanFruits")])
Iterations = 16005:26000
Thinning interval = 5
Number of chains = 3
Sample size per chain = 2000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
1060.4967 48.5300 0.6265 0.6341
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
969.4 1027.1 1058.3 1091.8 1160.8 Not too far from the true value which is 1000.
midpoint
plot(mcmc1[,c("midpoint")])
summary(mcmc1[,c("midpoint")])
Iterations = 16005:26000
Thinning interval = 5
Number of chains = 3
Sample size per chain = 2000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
529.3575 30.2216 0.3902 0.6382
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
471.5 503.9 535.6 557.9 565.1 Not too far from the true value of 500, but the shape of the distribution is a bit surprising…
steepness
plot(mcmc1[,c("steepness")])
summary(mcmc1[,c("steepness")])
Iterations = 16005:26000
Thinning interval = 5
Number of chains = 3
Sample size per chain = 2000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
1.89583 1.49891 0.01935 0.03273
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
0.08961 0.49946 1.57574 3.11968 4.78204 Far from the true value of 0.1! There may not be much information in the data to estimate the steepness parameter, or maybe there is a correlation between two parameters.
scalePrior
plot(mcmc1[,c("scalePrior")])
summary(mcmc1[,c("scalePrior")])
Iterations = 16005:26000
Thinning interval = 5
Number of chains = 3
Sample size per chain = 2000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
8.98131 1.66137 0.02145 0.02179
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
6.055 7.804 8.857 10.030 12.579 Not too far from the true value which was 10.
Are there correlations between parameters in our model?
crosscorr.plot(mcmc1)
There does not seem to be strong correlations between the variables.
More complex model that includes some weather variable
Relationship between precipitation at year n and the probability of making fruits at year n
We use a linear function of the absolute distance to the optimum amount of precipitation that returns values between 0 and 1.
abs_dist <- function(x, optimum=50) {
return (1-abs(x-50)/50)
}Demonstration of the link between precipitation and probability of making fruits
x=1:100
plot(abs_dist(x), t="l", ylab="Probability of making fruits", xlab="Precipitation")
Simulation that includes weather and influence of the previous year
# Simulates fructification at a given year with a model that depends on the previous year's number of fruits and the precipitations at year n.
# This dependency between years is modelled through the ref_logistic function, which provides a probability of making buds p_buds.
# Then we use abs_dist to get a probability p_1 to get fruits from buds, based on precipitation at year n.
simulatePrecipitations <- function(Nyears=7) {
return(runif(Nyears, min=0, max=100))
}
simulateFructification <- function (fruits_last_year, precipitation, midpoint, maximum=1, steepness=0.1, meanFruits=1000, theta_fruits=10) {
p_buds = ref_logistic(fruits_last_year, midpoint, maximum, steepness)
p1 = abs_dist(precipitation)
#print(paste("Probability of getting buds this year: ", p_buds, sep=""))
n_fruits = rnegbin(1, mu=meanFruits * p_buds * p1, theta=theta_fruits)
return(n_fruits)
}
# Simulates fructifications over N years, by calling iteratively simulateFructification.
simulateNYearsFructification <- function (Nyears=7, midpoint=500, maximum=1, steepness=0.1, meanFruits=1000, theta_fruits=10) {
# We simulate precipitations
precipitations = c(50)
precipitations = c(precipitations, simulatePrecipitations(Nyears))
fruits = c()
# we arbitrarily start assuming there had been fruits the year preceding the samples, with the best amount of precipitation
fruits = c(fruits, rnegbin(1, mu=meanFruits, theta=theta_fruits))
for (year in 1:Nyears) {
fruits = c(fruits, simulateFructification(fruits[length(fruits)], precipitations[year], midpoint, maximum, steepness, meanFruits, theta_fruits))
}
return(list(precipitations, fruits))
}
# Simulates fructifications over N years and in several areas, by calling iteratively simulateNYearsFructification.
simulateNAreas <- function(Nareas=15) {
all_fruits = c()
all_precipitations=c()
for (area in 1:Nareas) {
precipitations_fruits_li = simulateNYearsFructification()
all_precipitations = rbind(all_precipitations, precipitations_fruits_li[[1]])
all_fruits = rbind(all_fruits, precipitations_fruits_li[[2]])
}
return(list(all_precipitations, all_fruits))
}
Let’s look at one simulation
precipitations_fruits_li = simulateNYearsFructification()par(mfrow=c(2,1), mar=c(3, 4, 0, 2))
plot(x=1:8, precipitations_fruits_li[[1]], ylab="Precipitations", xlab="", pch=20)
plot(x=1:8, precipitations_fruits_li[[2]], ylab="Number of fruits", xlab="year", pch=20)
full simulation
precipitations_all_fruits_li <- simulateNAreas()
all_precipitations <- precipitations_all_fruits_li[[1]]
all_fruits <- precipitations_all_fruits_li[[2]]
Estimation under the complex model m2
jags_model_2 <-
"model
{
for (area in 1:Nareas)
{
# special case: year 1
fruits[area, 1] ~ dnegbin(min(scalePrior/(mu[area, 1]+scalePrior), 0.999), scalePrior)
mu[area, 1] <- pBuds[area, 1] * meanFruits * p1[area, 1]
pBuds[area, 1] ~ dbeta(1.0,1.0)
p1[area, 1] <- 1-abs(precipitation[area,1]-optimalPrecipitation)/(max(100-optimalPrecipitation, optimalPrecipitation))
# years 2 and following
for (year in 2:Nyears)
{
fruits[area, year] ~ dnegbin(min(scalePrior/(mu[area, year]+scalePrior), 0.999), scalePrior)
mu[area,year] <- pBuds[area, year] * meanFruits * p1[area, year]
pBuds[area, year] <- 1.0 - 1.0/(1+exp(-steepness*(fruits[area, year-1]-midpoint)))
p1[area, year] <- 1-abs(precipitation[area,1]-optimalPrecipitation)/(max(100-optimalPrecipitation, optimalPrecipitation))
}
}
optimalPrecipitation ~ dunif(0,100)
scalePrior ~ dunif(0.1, 5000)
meanFruits ~ dunif(5, 5000)
midpoint ~ dunif(0, 5000)
steepness ~ dunif(0, 5)
}"dcomplete = list(fruits=all_fruits, Nareas=dim(all_fruits)[1], Nyears=dim(all_fruits)[2], precipitation=all_precipitations)
ini <- list(list(optimalPrecipitation=10, scalePrior=1, meanFruits = 500, midpoint = 500, steepness = 0.1),
list(optimalPrecipitation=90, scalePrior=10, meanFruits = 1000, midpoint = 1000, steepness = 1.0),
list(optimalPrecipitation=40, scalePrior=100, meanFruits = 100, midpoint = 100, steepness = 4.0))
m2 <- jags.model(file=textConnection(jags_model_2), data = dcomplete, inits = ini, n.chains = 3)Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 120
Unobserved stochastic nodes: 20
Total graph size: 1026
Initializing model
|
| | 0%
|
|+ | 2%
|
|++ | 4%
|
|+++ | 6%
|
|++++ | 8%
|
|+++++ | 10%
|
|++++++ | 12%
|
|+++++++ | 14%
|
|++++++++ | 16%
|
|+++++++++ | 18%
|
|++++++++++ | 20%
|
|+++++++++++ | 22%
|
|++++++++++++ | 24%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 28%
|
|+++++++++++++++ | 30%
|
|++++++++++++++++ | 32%
|
|+++++++++++++++++ | 34%
|
|++++++++++++++++++ | 36%
|
|+++++++++++++++++++ | 38%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 42%
|
|++++++++++++++++++++++ | 44%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 48%
|
|+++++++++++++++++++++++++ | 50%
|
|++++++++++++++++++++++++++ | 52%
|
|+++++++++++++++++++++++++++ | 54%
|
|++++++++++++++++++++++++++++ | 56%
|
|+++++++++++++++++++++++++++++ | 58%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 62%
|
|++++++++++++++++++++++++++++++++ | 64%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 68%
|
|+++++++++++++++++++++++++++++++++++ | 70%
|
|++++++++++++++++++++++++++++++++++++ | 72%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|++++++++++++++++++++++++++++++++++++++ | 76%
|
|+++++++++++++++++++++++++++++++++++++++ | 78%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 82%
|
|++++++++++++++++++++++++++++++++++++++++++ | 84%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 88%
|
|+++++++++++++++++++++++++++++++++++++++++++++ | 90%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 92%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|++++++++++++++++++++++++++++++++++++++++++++++++ | 96%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 98%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
# burnin
update(m2, 5000)
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
# sampling
mcmc2 <- coda.samples(m2, c("optimalPrecipitation", "scalePrior", "meanFruits", "midpoint", "steepness"), n.iter = 10000, thin = 5)
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%Convergence diagnostics
gelman.diag(mcmc2)Potential scale reduction factors:
Point est. Upper C.I.
meanFruits 1.00 1.00
midpoint 1.01 1.02
optimalPrecipitation 1.03 1.09
scalePrior 1.00 1.00
steepness 1.14 1.27
Multivariate psrf
1.03Some potential scale reduction factors seem a bit far from 1.0. The steepness in particular may be a bit problematic to estimate.
How many independent iterations did we get?
effectiveSize(mcmc2) meanFruits midpoint optimalPrecipitation scalePrior steepness
1886.3595 199.2191 478.6848 5331.7645 131.4734 Clearly, midpoint and steepness are difficult to sample. More samples would be needed.
How many iterations would be needed?
raftery.diag(mcmc2) [[1]]
Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95
You need a sample size of at least 3746 with these values of q, r and s
[[2]]
Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95
You need a sample size of at least 3746 with these values of q, r and s
[[3]]
Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95
You need a sample size of at least 3746 with these values of q, r and sSo we don’t have enough iterations:
dim(mcmc2[[1]])[1][1] 2000mcmc2 <- coda.samples(m2, c("optimalPrecipitation", "scalePrior", "meanFruits", "midpoint", "steepness"), n.iter = 30000, thin = 5)
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%How are we doing now?
gelman.diag(mcmc2)Potential scale reduction factors:
Point est. Upper C.I.
meanFruits 1.00 1.00
midpoint 1.00 1.01
optimalPrecipitation 1.00 1.00
scalePrior 1.00 1.00
steepness 1.01 1.02
Multivariate psrf
1This diagnostic seems good.
raftery.diag(mcmc2)
This one too.
effectiveSize(mcmc2) meanFruits midpoint optimalPrecipitation scalePrior steepness
20070.6589 985.6415 4469.0865 25149.2075 465.1268 Midpoint and steepness are really difficult to estimate: even though we have sampled for 100 000 iterations, in the end the autocorrelation is so strong it’s as if we only had 465 independent draws from the posterior for steepness! To improve the estimation of these parameters, one could run the chains longer, or change the model, structurally, or with the priors.
Are there correlations between parameters in our model?
crosscorr.plot(mcmc2)
There is a correlation between two variables:
crosscorr(mcmc2) meanFruits midpoint optimalPrecipitation scalePrior steepness
meanFruits 1.00000000 -0.043662376 -0.012425025 -0.06231696 -0.013420287
midpoint -0.04366238 1.000000000 -0.002899814 -0.15313660 0.665963409
optimalPrecipitation -0.01242503 -0.002899814 1.000000000 0.01274707 -0.008446693
scalePrior -0.06231696 -0.153136600 0.012747075 1.00000000 -0.111870177
steepness -0.01342029 0.665963409 -0.008446693 -0.11187018 1.0000000000.67: no wonder steepness and midpoint are difficult to estimate: they are strongly correlated.
Analysis of the estimated parameters
plot(mcmc2)
We can see on the traces that steepness and midpoint are difficult to sample for the mcmc.
Estimated parameters
summary(mcmc2)
Iterations = 46005:146000
Thinning interval = 5
Number of chains = 3
Sample size per chain = 20000
1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:
Mean SD Naive SE Time-series SE
meanFruits 1137.4237 243.0871 0.9923990 1.715768
midpoint 517.7311 27.5409 0.1124353 0.994977
optimalPrecipitation 49.4480 31.2113 0.1274198 0.467075
scalePrior 1.2794 0.1949 0.0007956 0.001243
steepness 0.4591 0.9029 0.0036862 0.080577
2. Quantiles for each variable:
2.5% 25% 50% 75% 97.5%
meanFruits 711.68826 949.88293 1136.0055 1310.4696 1613.231
midpoint 476.04084 497.99014 511.5795 531.6963 572.150
optimalPrecipitation 1.83924 20.32795 48.8303 78.5571 98.117
scalePrior 0.93085 1.14210 1.2685 1.4036 1.694
steepness 0.06404 0.09414 0.1232 0.2003 3.773Looking at the median values, the estimates are pretty good for most parameters. The true values were 1000 for meanFruits, 500 for midpoint, 50 for optimalPrecipitation, 0.1 for steepness. But it was 10 for scalePrior, which is estimated here at 1.27. One would need to investigate why scalePrior is not well estimated.
LS0tCnRpdGxlOiAiQmFzaWMgbWFzdGluZyBtb2RlbCIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3J9CmxpYnJhcnkoTUFTUykKbGlicmFyeShyamFncykKYGBgCgoKIyBTaW11bGF0aW9uCgojIyBSZWxhdGlvbnNoaXAgYmV0d2VlbiB0aGUgbnVtYmVyIG9mIGZydWl0cyBhdCB5ZWFyIG4tMSBhbmQgdGhlIHByb2JhYmlsaXR5IG9mIG1ha2luZyBidWRzIGF0IHllYXIgbgpXZSB1c2UgYSBsb2dpc3RpYyBmdW5jdGlvbiBhbmQgcmVmbGVjdCBpdCB0byBtYWtlIGl0IGRlY3JlYXNlIHdpdGggaW5jcmVhc2luZyBmcnVpdHMgYXQgeWVhciBuLTEuCgpgYGB7cn0KIyBMb2dpc3RpYyBmdW5jdGlvbjogeSBnb2VzIGZyb20gMCB0byBtYXhpbXVtIGFzIHggaW5jcmVhc2VzCmxvZ2lzdGljIDwtIGZ1bmN0aW9uKHgsIG1pZHBvaW50LCBtYXhpbXVtPTEsIHN0ZWVwbmVzcz0wLjEpIHsKICB5ID0gbWF4aW11bS8oMStleHAoLXN0ZWVwbmVzcyooeC1taWRwb2ludCkpKQogIHJldHVybih5KQp9CgojIFJlZmxlY3RlZCBsb2dpc3RpYyBmdW5jdGlvbjogeSBnb2VzIGZyb20gbWF4aW11bSB0byAwIGFzIHggaW5jcmVhc2VzCnJlZl9sb2dpc3RpYyA8LSBmdW5jdGlvbih4LCBtaWRwb2ludCwgbWF4aW11bT0xLCBzdGVlcG5lc3M9MC4xKSB7CiByZXR1cm4obWF4aW11bSAtIGxvZ2lzdGljKHgsIG1pZHBvaW50LCBtYXhpbXVtLCBzdGVlcG5lc3MpKSAKfQoKYGBgCgoKIyMjIFBsb3R0aW5nIHRoZSBsb2dpc3RpYyBmdW5jdGlvbiB3ZSB1c2UKYGBge3J9CnggPSAoMToxMDAwKSoyCnBsb3QoeCwgcmVmX2xvZ2lzdGljKHgsIG1pZHBvaW50PTUwMCksIHQ9ImwiKSAgCmBgYAoKCmBgYHtyfQojIFNpbXVsYXRlcyBmcnVjdGlmaWNhdGlvbiBhdCBhIGdpdmVuIHllYXIgd2l0aCBhIHNpbXBsZSBtb2RlbCB0aGF0IGRlcGVuZHMgb24gdGhlIHByZXZpb3VzIHllYXIncyBudW1iZXIgb2YgZnJ1aXRzIG9ubHkuIAojIFRoaXMgZGVwZW5kZW5jeSBpcyBtb2RlbGxlZCB0aHJvdWdoIHRoZSByZWZfbG9naXN0aWMgZnVuY3Rpb24sIHdoaWNoIHByb3ZpZGVzIGEgcHJvYmFiaWxpdHkgb2YgbWFraW5nIGJ1ZHMgcF9idWRzLgpzaW11bGF0ZUZydWN0aWZpY2F0aW9uIDwtIGZ1bmN0aW9uIChmcnVpdHNfbGFzdF95ZWFyLCBtaWRwb2ludCwgbWF4aW11bT0xLCBzdGVlcG5lc3M9MC4xLCBtZWFuRnJ1aXRzPTEwMDAsIHRoZXRhX2ZydWl0cz0xMCkgewogIHBfYnVkcyA9IHJlZl9sb2dpc3RpYyhmcnVpdHNfbGFzdF95ZWFyLCBtaWRwb2ludCwgbWF4aW11bSwgc3RlZXBuZXNzKQogICNwcmludChwYXN0ZSgiUHJvYmFiaWxpdHkgb2YgZ2V0dGluZyBidWRzIHRoaXMgeWVhcjogIiwgcF9idWRzLCBzZXA9IiIpKQogIG5fZnJ1aXRzID0gcm5lZ2JpbigxLCBtdT1tZWFuRnJ1aXRzICogcF9idWRzLCB0aGV0YT10aGV0YV9mcnVpdHMpCiAgcmV0dXJuKG5fZnJ1aXRzKQp9CgojIFNpbXVsYXRlcyBmcnVjdGlmaWNhdGlvbnMgb3ZlciBOIHllYXJzLCBieSBjYWxsaW5nIGl0ZXJhdGl2ZWx5IHNpbXVsYXRlRnJ1Y3RpZmljYXRpb24uCnNpbXVsYXRlTlllYXJzRnJ1Y3RpZmljYXRpb24gPC0gZnVuY3Rpb24gKE55ZWFycz03LCBtaWRwb2ludD01MDAsIG1heGltdW09MSwgc3RlZXBuZXNzPTAuMSwgbWVhbkZydWl0cz0xMDAwLCB0aGV0YV9mcnVpdHM9MTApIHsKICBmcnVpdHMgPSBjKCkKICAjIHdlIGFyYml0cmFyaWx5IHN0YXJ0IGFzc3VtaW5nIHRoZXJlIGhhZCBiZWVuIGZydWl0cyB0aGUgeWVhciBwcmVjZWRpbmcgdGhlIHNhbXBsZXMKICBmcnVpdHMgPSBjKGZydWl0cywgcm5lZ2JpbigxLCBtdT1tZWFuRnJ1aXRzLCB0aGV0YT10aGV0YV9mcnVpdHMpKQogICAgZm9yICh5ZWFyIGluIDE6TnllYXJzKSB7CiAgICAgIGZydWl0cyA9IGMoZnJ1aXRzLCBzaW11bGF0ZUZydWN0aWZpY2F0aW9uKGZydWl0c1tsZW5ndGgoZnJ1aXRzKV0sIG1pZHBvaW50LCBtYXhpbXVtLCBzdGVlcG5lc3MsIG1lYW5GcnVpdHMsIHRoZXRhX2ZydWl0cykpCiAgICB9CiAgcmV0dXJuKGZydWl0cykKfQoKIyBTaW11bGF0ZXMgZnJ1Y3RpZmljYXRpb25zIG92ZXIgTiB5ZWFycyBhbmQgaW4gc2V2ZXJhbCBhcmVhcywgYnkgY2FsbGluZyBpdGVyYXRpdmVseSBzaW11bGF0ZU5ZZWFyc0ZydWN0aWZpY2F0aW9uLgpzaW11bGF0ZU5BcmVhcyA8LSBmdW5jdGlvbihOYXJlYXM9MTUpIHsKICBhbGxfZnJ1aXRzID0gYygpCiAgZm9yIChhcmVhIGluIDE6TmFyZWFzKSB7CiAgICBhbGxfZnJ1aXRzID0gcmJpbmQoYWxsX2ZydWl0cywgc2ltdWxhdGVOWWVhcnNGcnVjdGlmaWNhdGlvbigpKSAgCiAgfQogIHJldHVybihhbGxfZnJ1aXRzKQp9CgpgYGAKCiMjIFNpbXVsYXRpbmcgb25lIGFyZWEKYGBge3J9CmZydWl0cyA9IHNpbXVsYXRlTlllYXJzRnJ1Y3RpZmljYXRpb24oKQpgYGAKCiMjIFdoYXQgZG9lcyBpdCBsb29rIGxpa2U/CgpgYGB7cn0KcGxvdChmcnVpdHMsIHhsYWI9IlllYXIiKQpgYGAKCiMjIFNpbXVsYXRpbmcgMTUgYXJlYXMKYGBge3J9CmZydWl0cyA9IHNpbXVsYXRlTkFyZWFzKCkKYGBgCgoKIyMgV2hhdCBkb2VzIGl0IGxvb2sgbGlrZT8KCmBgYHtyfQpwbG90KHg9MTo4LCB5PWZydWl0c1sxLF0sIHhsYWI9IlllYXIiLCB5bGFiPSJOdW1iZXIgb2YgZnJ1aXRzIiwgdD0ibCIsIGNvbD0wLCB5bGltPWMoMCwyMDAwKSkKZm9yIChpIGluIDI6ZGltKGZydWl0cylbMV0pIHsKICBsaW5lcyhmcnVpdHNbaSxdLCBjb2w9aSkKfQpgYGAKCiMgRXN0aW1hdGlvbgoKV2UgYXJlIGdvaW5nIHRvIHVzZSBKYWdzIHRvIGVzdGltYXRlIHRoZSBwYXJhbWV0ZXJzIG9mIGEgbW9kZWwgdGhhdCBtYXRjaGVzIHRoZSBzaW11bGF0aW9uIG1vZGVsLgpUaGUgZGlzdHJpYnV0aW9uIGRuZWdiaW4gaW4gSmFncyBpcyBwYXJhbWV0cml6ZWQgaW4gcCBhbmQgcywgd2hlcmVhcyBSJ3Mgcm5lZ2JpbiBpcyBwYXJhbWV0cml6ZWQgaW4gbXUgKHRoZSBtZWFuKSBhbmQgdGhldGEuClRoZSByZWxhdGlvbnNoaXBzIGJldHdlZW4gdGhvc2UgcGFyYW1ldGVycyBhcmUgYXMgZm9sbG93czoKIC0gIHMgPSB0aGV0YQogLSAgbXUgPSBzKDEtcCkvcCBpLmUuIHAgPSBzLyhtdStzKQogSW4gdGhlIHNpbXVsYXRpb24sIHdlIGhhdmU6IG11ID0gcF9idWRzICogcDEgKiBtZWFuRnJ1aXRzCldlIHdhbnQgdG8gaGF2ZSB0aGVzZSBwYXJhbWV0ZXJzIGluIHRoZSBKYWdzIGVzdGltYXRpb24gbW9kZWwgdG8gYmUgYWJsZSB0byBjb21wYXJlIHRoZSBlc3RpbWF0ZWQgdG8gdGhlIHRydWUgdmFsdWVzOiB0aGUgZXN0aW1hdGlvbiBtb2RlbCBoYXMgdG8gbWF0Y2ggdGhlIHNpbXVsYXRpb24gbW9kZWwuCiAKYGBge3J9CgpqYWdzX21vZGVsXzEgPC0KIm1vZGVsCnsKZm9yIChhcmVhIGluIDE6TmFyZWFzKQp7CiAgIyBzcGVjaWFsIGNhc2U6IHllYXIgMQogIGZydWl0c1thcmVhLCAxXSB+IGRuZWdiaW4obWluKHNjYWxlUHJpb3IvKG11W2FyZWEsIDFdK3NjYWxlUHJpb3IpLCAwLjk5OSksIHNjYWxlUHJpb3IpCiAgbXVbYXJlYSwgMV0gPC0gcEJ1ZHNbYXJlYSwgMV0gKiBtZWFuRnJ1aXRzCiAgcEJ1ZHNbYXJlYSwgMV0gfiBkYmV0YSgxLjAsMS4wKQogICMgeWVhcnMgMiBhbmQgZm9sbG93aW5nCiAgZm9yICh5ZWFyIGluIDI6TnllYXJzKQogIHsKICAgIGZydWl0c1thcmVhLCB5ZWFyXSB+IGRuZWdiaW4obWluKHNjYWxlUHJpb3IvKG11W2FyZWEsIHllYXJdK3NjYWxlUHJpb3IpLCAwLjk5OSksICBzY2FsZVByaW9yKQogICAgbXVbYXJlYSx5ZWFyXSA8LSBwQnVkc1thcmVhLCB5ZWFyXSAqIG1lYW5GcnVpdHMKICAgIHBCdWRzW2FyZWEsIHllYXJdIDwtIDEuMCAtIDEuMC8oMStleHAoLXN0ZWVwbmVzcyooZnJ1aXRzW2FyZWEsIHllYXItMV0tbWlkcG9pbnQpKSkKICB9Cn0KICBzY2FsZVByaW9yIH4gZHVuaWYoMC4xLCA1MDAwKQogIG1lYW5GcnVpdHMgfiBkdW5pZig1LCA1MDAwKQogIG1pZHBvaW50IH4gZHVuaWYoMCwgNTAwMCkKICBzdGVlcG5lc3MgfiBkdW5pZigwLCA1KQp9IgpgYGAKCgpgYGB7cn0KZGNvbXBsZXRlID0gbGlzdChmcnVpdHM9ZnJ1aXRzLCBOYXJlYXM9ZGltKGZydWl0cylbMV0sIE55ZWFycz1kaW0oZnJ1aXRzKVsyXSkKCmluaSA8LSBsaXN0KGxpc3Qoc2NhbGVQcmlvcj0xLCBtZWFuRnJ1aXRzID0gNTAwLCBtaWRwb2ludCA9IDUwMCwgc3RlZXBuZXNzID0gMC4xKSwKbGlzdChzY2FsZVByaW9yPTEwLCBtZWFuRnJ1aXRzID0gMTAwMCwgbWlkcG9pbnQgPSAxMDAwLCBzdGVlcG5lc3MgPSAxLjApLApsaXN0KHNjYWxlUHJpb3I9MTAwLCBtZWFuRnJ1aXRzID0gMTAwLCBtaWRwb2ludCA9IDEwMCwgc3RlZXBuZXNzID0gNC4wKSkKbTE8LWphZ3MubW9kZWwoZmlsZSA9IHRleHRDb25uZWN0aW9uKGphZ3NfbW9kZWxfMSksIGRhdGEgPSBkY29tcGxldGUsIGluaXRzID0gaW5pLCBuLmNoYWlucyA9IDMpCgpgYGAKCgpgYGB7cn0KCiMgYnVybmluCnVwZGF0ZShtMSwgNTAwMCkKCmBgYAoKCmBgYHtyfQoKIyBzYW1wbGluZwptY21jMSA8LSBjb2RhLnNhbXBsZXMobTEsIGMoIm1lYW5GcnVpdHMiLCAibWlkcG9pbnQiLCAic3RlZXBuZXNzIiwgInNjYWxlUHJpb3IiKSwgbi5pdGVyID0gMTAwMDAsIHRoaW4gPSA1KQoKYGBgCgoKIyMgQ29udmVyZ2VuY2UgZGlhZ25vc3RpY3MKCmBgYHtyfQpnZWxtYW4uZGlhZyhtY21jMSkKYGBgCgpUaGUgcG90ZW50aWFsIHNjYWxlIHJlZHVjdGlvbiBmYWN0b3JzIHNlZW0gY2xvc2UgdG8gMS4wIGZvciBib3RoIHRoZSBtZWFuIGFuZCB0aGUgdXBwZXIgQy5JLiwgd2hpY2ggaXMgYSBnb29kIHNpZ24gdGhhdCBhbGwgMyBjaGFpbnMgY29udmVyZ2VkIHRvIHNpbWlsYXIgcG9zdGVyaW9yIGRlbnNpdGllcy4KCmBgYHtyfQplZmZlY3RpdmVTaXplKG1jbWMxKQpgYGAKCkl0IHdvdWxkIHNlZW0gbGlrZSB3ZSBoYXZlIGVub3VnaCBpbmRlcGVuZGVudCBzYW1wbGVzIHRvIGVzdGltYXRlIHRoZSBwb3N0ZXJpb3IgZGlzdHJpYnV0aW9ucyB3ZWxsLgoKCiMjIEFuYWx5c2lzIG9mIHRoZSBlc3RpbWF0ZWQgcGFyYW1ldGVycwoKYGBge3J9CgpwbG90KG1jbWMxKQpgYGAKCldlIGRvbid0IHNlZSBtdWNoLCB3ZSdsbCBsb29rIGF0IGVhY2ggdmFyaWFibGUgaW4gdHVybi4KCiMjIyBtZWFuRnJ1aXRzCmBgYHtyfQpwbG90KG1jbWMxWyxjKCJtZWFuRnJ1aXRzIildKQpgYGAKCgpgYGB7cn0Kc3VtbWFyeShtY21jMVssYygibWVhbkZydWl0cyIpXSkKYGBgCgpOb3QgdG9vIGZhciBmcm9tIHRoZSB0cnVlIHZhbHVlIHdoaWNoIGlzIDEwMDAuCgojIyMgbWlkcG9pbnQKCmBgYHtyfQpwbG90KG1jbWMxWyxjKCJtaWRwb2ludCIpXSkKYGBgCgpgYGB7cn0Kc3VtbWFyeShtY21jMVssYygibWlkcG9pbnQiKV0pCmBgYAoKTm90IHRvbyBmYXIgZnJvbSB0aGUgdHJ1ZSB2YWx1ZSBvZiA1MDAsIGJ1dCB0aGUgc2hhcGUgb2YgdGhlIGRpc3RyaWJ1dGlvbiBpcyBhIGJpdCBzdXJwcmlzaW5nLi4uCgoKIyMjIHN0ZWVwbmVzcwoKYGBge3J9CnBsb3QobWNtYzFbLGMoInN0ZWVwbmVzcyIpXSkKYGBgCgoKYGBge3J9CnN1bW1hcnkobWNtYzFbLGMoInN0ZWVwbmVzcyIpXSkKYGBgCgpGYXIgZnJvbSB0aGUgdHJ1ZSB2YWx1ZSBvZiAwLjEhIFRoZXJlIG1heSBub3QgYmUgbXVjaCBpbmZvcm1hdGlvbiBpbiB0aGUgZGF0YSB0byBlc3RpbWF0ZSB0aGUgc3RlZXBuZXNzIHBhcmFtZXRlciwgb3IgbWF5YmUgdGhlcmUgaXMgYSBjb3JyZWxhdGlvbiBiZXR3ZWVuIHR3byBwYXJhbWV0ZXJzLgoKIyMjIHNjYWxlUHJpb3IKYGBge3J9CnBsb3QobWNtYzFbLGMoInNjYWxlUHJpb3IiKV0pCmBgYAoKCmBgYHtyfQpzdW1tYXJ5KG1jbWMxWyxjKCJzY2FsZVByaW9yIildKQpgYGAKCk5vdCB0b28gZmFyIGZyb20gdGhlIHRydWUgdmFsdWUgd2hpY2ggd2FzIDEwLgoKIyMgQXJlIHRoZXJlIGNvcnJlbGF0aW9ucyBiZXR3ZWVuIHBhcmFtZXRlcnMgaW4gb3VyIG1vZGVsPwpgYGB7cn0KY3Jvc3Njb3JyLnBsb3QobWNtYzEpCmBgYApUaGVyZSBkb2VzIG5vdCBzZWVtIHRvIGJlIHN0cm9uZyBjb3JyZWxhdGlvbnMgYmV0d2VlbiB0aGUgdmFyaWFibGVzLiAKCiMgTW9yZSBjb21wbGV4IG1vZGVsIHRoYXQgaW5jbHVkZXMgc29tZSB3ZWF0aGVyIHZhcmlhYmxlCgojIyBSZWxhdGlvbnNoaXAgYmV0d2VlbiBwcmVjaXBpdGF0aW9uIGF0IHllYXIgbiBhbmQgdGhlIHByb2JhYmlsaXR5IG9mIG1ha2luZyBmcnVpdHMgYXQgeWVhciBuCldlIHVzZSBhIGxpbmVhciBmdW5jdGlvbiBvZiB0aGUgYWJzb2x1dGUgZGlzdGFuY2UgdG8gdGhlIG9wdGltdW0gYW1vdW50IG9mIHByZWNpcGl0YXRpb24gdGhhdCByZXR1cm5zIHZhbHVlcyBiZXR3ZWVuIDAgYW5kIDEuIAoKYGBge3J9CmFic19kaXN0IDwtIGZ1bmN0aW9uKHgsIG9wdGltdW09NTApIHsKICByZXR1cm4gKDEtYWJzKHgtNTApLzUwKQp9CmBgYAoKIyMjIERlbW9uc3RyYXRpb24gb2YgdGhlIGxpbmsgYmV0d2VlbiBwcmVjaXBpdGF0aW9uIGFuZCBwcm9iYWJpbGl0eSBvZiBtYWtpbmcgZnJ1aXRzCmBgYHtyfQp4PTE6MTAwCnBsb3QoYWJzX2Rpc3QoeCksIHQ9ImwiLCB5bGFiPSJQcm9iYWJpbGl0eSBvZiBtYWtpbmcgZnJ1aXRzIiwgeGxhYj0iUHJlY2lwaXRhdGlvbiIpCmBgYAoKIyMgU2ltdWxhdGlvbiB0aGF0IGluY2x1ZGVzIHdlYXRoZXIgYW5kIGluZmx1ZW5jZSBvZiB0aGUgcHJldmlvdXMgeWVhcgoKYGBge3J9CiMgU2ltdWxhdGVzIGZydWN0aWZpY2F0aW9uIGF0IGEgZ2l2ZW4geWVhciB3aXRoIGEgbW9kZWwgdGhhdCBkZXBlbmRzIG9uIHRoZSBwcmV2aW91cyB5ZWFyJ3MgbnVtYmVyIG9mIGZydWl0cyBhbmQgdGhlIHByZWNpcGl0YXRpb25zIGF0IHllYXIgbi4gCiMgVGhpcyBkZXBlbmRlbmN5IGJldHdlZW4geWVhcnMgaXMgbW9kZWxsZWQgdGhyb3VnaCB0aGUgcmVmX2xvZ2lzdGljIGZ1bmN0aW9uLCB3aGljaCBwcm92aWRlcyBhIHByb2JhYmlsaXR5IG9mIG1ha2luZyBidWRzIHBfYnVkcy4KIyBUaGVuIHdlIHVzZSBhYnNfZGlzdCB0byBnZXQgYSBwcm9iYWJpbGl0eSBwXzEgdG8gZ2V0IGZydWl0cyBmcm9tIGJ1ZHMsIGJhc2VkIG9uIHByZWNpcGl0YXRpb24gYXQgeWVhciBuLgoKc2ltdWxhdGVQcmVjaXBpdGF0aW9ucyA8LSBmdW5jdGlvbihOeWVhcnM9NykgewogIHJldHVybihydW5pZihOeWVhcnMsIG1pbj0wLCBtYXg9MTAwKSkKfQoKc2ltdWxhdGVGcnVjdGlmaWNhdGlvbiA8LSBmdW5jdGlvbiAoZnJ1aXRzX2xhc3RfeWVhciwgcHJlY2lwaXRhdGlvbiwgbWlkcG9pbnQsIG1heGltdW09MSwgc3RlZXBuZXNzPTAuMSwgbWVhbkZydWl0cz0xMDAwLCB0aGV0YV9mcnVpdHM9MTApIHsKICBwX2J1ZHMgPSByZWZfbG9naXN0aWMoZnJ1aXRzX2xhc3RfeWVhciwgbWlkcG9pbnQsIG1heGltdW0sIHN0ZWVwbmVzcykKICBwMSA9IGFic19kaXN0KHByZWNpcGl0YXRpb24pCiAgI3ByaW50KHBhc3RlKCJQcm9iYWJpbGl0eSBvZiBnZXR0aW5nIGJ1ZHMgdGhpcyB5ZWFyOiAiLCBwX2J1ZHMsIHNlcD0iIikpCiAgbl9mcnVpdHMgPSBybmVnYmluKDEsIG11PW1lYW5GcnVpdHMgKiBwX2J1ZHMgKiBwMSwgdGhldGE9dGhldGFfZnJ1aXRzKQogIHJldHVybihuX2ZydWl0cykKfQoKIyBTaW11bGF0ZXMgZnJ1Y3RpZmljYXRpb25zIG92ZXIgTiB5ZWFycywgYnkgY2FsbGluZyBpdGVyYXRpdmVseSBzaW11bGF0ZUZydWN0aWZpY2F0aW9uLgpzaW11bGF0ZU5ZZWFyc0ZydWN0aWZpY2F0aW9uIDwtIGZ1bmN0aW9uIChOeWVhcnM9NywgbWlkcG9pbnQ9NTAwLCBtYXhpbXVtPTEsIHN0ZWVwbmVzcz0wLjEsIG1lYW5GcnVpdHM9MTAwMCwgdGhldGFfZnJ1aXRzPTEwKSB7CiAgIyBXZSBzaW11bGF0ZSBwcmVjaXBpdGF0aW9ucwogIHByZWNpcGl0YXRpb25zID0gYyg1MCkKICBwcmVjaXBpdGF0aW9ucyA9IGMocHJlY2lwaXRhdGlvbnMsIHNpbXVsYXRlUHJlY2lwaXRhdGlvbnMoTnllYXJzKSkKICBmcnVpdHMgPSBjKCkKICAjIHdlIGFyYml0cmFyaWx5IHN0YXJ0IGFzc3VtaW5nIHRoZXJlIGhhZCBiZWVuIGZydWl0cyB0aGUgeWVhciBwcmVjZWRpbmcgdGhlIHNhbXBsZXMsIHdpdGggdGhlIGJlc3QgYW1vdW50IG9mIHByZWNpcGl0YXRpb24KICBmcnVpdHMgPSBjKGZydWl0cywgcm5lZ2JpbigxLCBtdT1tZWFuRnJ1aXRzLCB0aGV0YT10aGV0YV9mcnVpdHMpKQogICAgZm9yICh5ZWFyIGluIDE6TnllYXJzKSB7CiAgICAgIGZydWl0cyA9IGMoZnJ1aXRzLCBzaW11bGF0ZUZydWN0aWZpY2F0aW9uKGZydWl0c1tsZW5ndGgoZnJ1aXRzKV0sIHByZWNpcGl0YXRpb25zW3llYXJdLCBtaWRwb2ludCwgbWF4aW11bSwgc3RlZXBuZXNzLCBtZWFuRnJ1aXRzLCB0aGV0YV9mcnVpdHMpKQogICAgfQogIHJldHVybihsaXN0KHByZWNpcGl0YXRpb25zLCBmcnVpdHMpKQp9CgojIFNpbXVsYXRlcyBmcnVjdGlmaWNhdGlvbnMgb3ZlciBOIHllYXJzIGFuZCBpbiBzZXZlcmFsIGFyZWFzLCBieSBjYWxsaW5nIGl0ZXJhdGl2ZWx5IHNpbXVsYXRlTlllYXJzRnJ1Y3RpZmljYXRpb24uCnNpbXVsYXRlTkFyZWFzIDwtIGZ1bmN0aW9uKE5hcmVhcz0xNSkgewogIGFsbF9mcnVpdHMgPSBjKCkKICBhbGxfcHJlY2lwaXRhdGlvbnM9YygpCiAgZm9yIChhcmVhIGluIDE6TmFyZWFzKSB7CiAgICBwcmVjaXBpdGF0aW9uc19mcnVpdHNfbGkgPSBzaW11bGF0ZU5ZZWFyc0ZydWN0aWZpY2F0aW9uKCkKICAgIGFsbF9wcmVjaXBpdGF0aW9ucyA9IHJiaW5kKGFsbF9wcmVjaXBpdGF0aW9ucywgcHJlY2lwaXRhdGlvbnNfZnJ1aXRzX2xpW1sxXV0pCiAgICBhbGxfZnJ1aXRzID0gcmJpbmQoYWxsX2ZydWl0cywgcHJlY2lwaXRhdGlvbnNfZnJ1aXRzX2xpW1syXV0pCiAgfQogIHJldHVybihsaXN0KGFsbF9wcmVjaXBpdGF0aW9ucywgYWxsX2ZydWl0cykpCn0KCgpgYGAKCiMjIExldCdzIGxvb2sgYXQgb25lIHNpbXVsYXRpb24KYGBge3J9CnByZWNpcGl0YXRpb25zX2ZydWl0c19saSA9IHNpbXVsYXRlTlllYXJzRnJ1Y3RpZmljYXRpb24oKQpgYGAKCmBgYHtyfQpwYXIobWZyb3c9YygyLDEpLCBtYXI9YygzLCA0LCAwLCAyKSkKcGxvdCh4PTE6OCwgcHJlY2lwaXRhdGlvbnNfZnJ1aXRzX2xpW1sxXV0sIHlsYWI9IlByZWNpcGl0YXRpb25zIiwgeGxhYj0iIiwgcGNoPTIwKQpwbG90KHg9MTo4LCBwcmVjaXBpdGF0aW9uc19mcnVpdHNfbGlbWzJdXSwgeWxhYj0iTnVtYmVyIG9mIGZydWl0cyIsIHhsYWI9InllYXIiLCBwY2g9MjApCgpgYGAKCiMjIGZ1bGwgc2ltdWxhdGlvbgpgYGB7cn0KcHJlY2lwaXRhdGlvbnNfYWxsX2ZydWl0c19saSA8LSBzaW11bGF0ZU5BcmVhcygpCmFsbF9wcmVjaXBpdGF0aW9ucyA8LSBwcmVjaXBpdGF0aW9uc19hbGxfZnJ1aXRzX2xpW1sxXV0KYWxsX2ZydWl0cyA8LSBwcmVjaXBpdGF0aW9uc19hbGxfZnJ1aXRzX2xpW1syXV0KCmBgYAoKCiMgRXN0aW1hdGlvbiB1bmRlciB0aGUgY29tcGxleCBtb2RlbCBtMgoKIApgYGB7cn0KCmphZ3NfbW9kZWxfMiA8LQoibW9kZWwKewpmb3IgKGFyZWEgaW4gMTpOYXJlYXMpCnsKICAjIHNwZWNpYWwgY2FzZTogeWVhciAxCiAgZnJ1aXRzW2FyZWEsIDFdIH4gZG5lZ2JpbihtaW4oc2NhbGVQcmlvci8obXVbYXJlYSwgMV0rc2NhbGVQcmlvciksIDAuOTk5KSwgc2NhbGVQcmlvcikKICBtdVthcmVhLCAxXSA8LSBwQnVkc1thcmVhLCAxXSAqIG1lYW5GcnVpdHMgKiBwMVthcmVhLCAxXQogIHBCdWRzW2FyZWEsIDFdIH4gZGJldGEoMS4wLDEuMCkKICBwMVthcmVhLCAxXSA8LSAxLWFicyhwcmVjaXBpdGF0aW9uW2FyZWEsMV0tb3B0aW1hbFByZWNpcGl0YXRpb24pLyhtYXgoMTAwLW9wdGltYWxQcmVjaXBpdGF0aW9uLCBvcHRpbWFsUHJlY2lwaXRhdGlvbikpCiAgIyB5ZWFycyAyIGFuZCBmb2xsb3dpbmcKICBmb3IgKHllYXIgaW4gMjpOeWVhcnMpCiAgewogICAgZnJ1aXRzW2FyZWEsIHllYXJdIH4gZG5lZ2JpbihtaW4oc2NhbGVQcmlvci8obXVbYXJlYSwgeWVhcl0rc2NhbGVQcmlvciksIDAuOTk5KSwgIHNjYWxlUHJpb3IpCiAgICBtdVthcmVhLHllYXJdIDwtIHBCdWRzW2FyZWEsIHllYXJdICogbWVhbkZydWl0cyAqIHAxW2FyZWEsIHllYXJdCiAgICBwQnVkc1thcmVhLCB5ZWFyXSA8LSAxLjAgLSAxLjAvKDErZXhwKC1zdGVlcG5lc3MqKGZydWl0c1thcmVhLCB5ZWFyLTFdLW1pZHBvaW50KSkpCiAgICBwMVthcmVhLCB5ZWFyXSA8LSAxLWFicyhwcmVjaXBpdGF0aW9uW2FyZWEsMV0tb3B0aW1hbFByZWNpcGl0YXRpb24pLyhtYXgoMTAwLW9wdGltYWxQcmVjaXBpdGF0aW9uLCBvcHRpbWFsUHJlY2lwaXRhdGlvbikpCiAgfQp9CiAgb3B0aW1hbFByZWNpcGl0YXRpb24gfiBkdW5pZigwLDEwMCkKICBzY2FsZVByaW9yIH4gZHVuaWYoMC4xLCA1MDAwKQogIG1lYW5GcnVpdHMgfiBkdW5pZig1LCA1MDAwKQogIG1pZHBvaW50IH4gZHVuaWYoMCwgNTAwMCkKICBzdGVlcG5lc3MgfiBkdW5pZigwLCA1KQp9IgpgYGAKCgpgYGB7cn0KZGNvbXBsZXRlID0gbGlzdChmcnVpdHM9YWxsX2ZydWl0cywgTmFyZWFzPWRpbShhbGxfZnJ1aXRzKVsxXSwgTnllYXJzPWRpbShhbGxfZnJ1aXRzKVsyXSwgcHJlY2lwaXRhdGlvbj1hbGxfcHJlY2lwaXRhdGlvbnMpCgppbmkgPC0gbGlzdChsaXN0KG9wdGltYWxQcmVjaXBpdGF0aW9uPTEwLCBzY2FsZVByaW9yPTEsIG1lYW5GcnVpdHMgPSA1MDAsIG1pZHBvaW50ID0gNTAwLCBzdGVlcG5lc3MgPSAwLjEpLApsaXN0KG9wdGltYWxQcmVjaXBpdGF0aW9uPTkwLCBzY2FsZVByaW9yPTEwLCBtZWFuRnJ1aXRzID0gMTAwMCwgbWlkcG9pbnQgPSAxMDAwLCBzdGVlcG5lc3MgPSAxLjApLApsaXN0KG9wdGltYWxQcmVjaXBpdGF0aW9uPTQwLCBzY2FsZVByaW9yPTEwMCwgbWVhbkZydWl0cyA9IDEwMCwgbWlkcG9pbnQgPSAxMDAsIHN0ZWVwbmVzcyA9IDQuMCkpCm0yIDwtIGphZ3MubW9kZWwoZmlsZT10ZXh0Q29ubmVjdGlvbihqYWdzX21vZGVsXzIpLCBkYXRhID0gZGNvbXBsZXRlLCBpbml0cyA9IGluaSwgbi5jaGFpbnMgPSAzKQoKYGBgCgoKYGBge3J9CgojIGJ1cm5pbgp1cGRhdGUobTIsIDUwMDApCgpgYGAKCgpgYGB7cn0KCiMgc2FtcGxpbmcKbWNtYzIgPC0gY29kYS5zYW1wbGVzKG0yLCBjKCJvcHRpbWFsUHJlY2lwaXRhdGlvbiIsICJzY2FsZVByaW9yIiwgIm1lYW5GcnVpdHMiLCAibWlkcG9pbnQiLCAic3RlZXBuZXNzIiksIG4uaXRlciA9IDEwMDAwLCB0aGluID0gNSkKCmBgYAoKCiMjIENvbnZlcmdlbmNlIGRpYWdub3N0aWNzCgpgYGB7cn0KZ2VsbWFuLmRpYWcobWNtYzIpCmBgYAoKU29tZSBwb3RlbnRpYWwgc2NhbGUgcmVkdWN0aW9uIGZhY3RvcnMgc2VlbSBhIGJpdCBmYXIgZnJvbSAxLjAuIFRoZSBzdGVlcG5lc3MgaW4gcGFydGljdWxhciBtYXkgYmUgYSBiaXQgcHJvYmxlbWF0aWMgdG8gZXN0aW1hdGUuCgojIyMgSG93IG1hbnkgaW5kZXBlbmRlbnQgaXRlcmF0aW9ucyBkaWQgd2UgZ2V0PwpgYGB7cn0KZWZmZWN0aXZlU2l6ZShtY21jMikKYGBgCgpDbGVhcmx5LCBtaWRwb2ludCBhbmQgc3RlZXBuZXNzIGFyZSBkaWZmaWN1bHQgdG8gc2FtcGxlLiBNb3JlIHNhbXBsZXMgd291bGQgYmUgbmVlZGVkLgoKIyMjIEhvdyBtYW55IGl0ZXJhdGlvbnMgd291bGQgYmUgbmVlZGVkPwpgYGB7cn0KcmFmdGVyeS5kaWFnKG1jbWMyKSAKCmBgYAoKU28gd2UgZG9uJ3QgaGF2ZSBlbm91Z2ggaXRlcmF0aW9uczoKYGBge3J9CmRpbShtY21jMltbMV1dKVsxXQpgYGAKCmBgYHtyfQptY21jMiA8LSBjb2RhLnNhbXBsZXMobTIsIGMoIm9wdGltYWxQcmVjaXBpdGF0aW9uIiwgInNjYWxlUHJpb3IiLCAibWVhbkZydWl0cyIsICJtaWRwb2ludCIsICJzdGVlcG5lc3MiKSwgbi5pdGVyID0gMTAwMDAwLCB0aGluID0gNSkKCmBgYAoKSG93IGFyZSB3ZSBkb2luZyBub3c/CgpgYGB7cn0KZ2VsbWFuLmRpYWcobWNtYzIpCmBgYAoKVGhpcyBkaWFnbm9zdGljIHNlZW1zIGdvb2QuCgpgYGB7cn0KcmFmdGVyeS5kaWFnKG1jbWMyKSAKCmBgYAoKVGhpcyBvbmUgdG9vLgoKCmBgYHtyfQplZmZlY3RpdmVTaXplKG1jbWMyKQpgYGAKCk1pZHBvaW50IGFuZCBzdGVlcG5lc3MgYXJlIHJlYWxseSBkaWZmaWN1bHQgdG8gZXN0aW1hdGU6IGV2ZW4gdGhvdWdoIHdlIGhhdmUgc2FtcGxlZCBmb3IgMTAwIDAwMCBpdGVyYXRpb25zLCBpbiB0aGUgZW5kIHRoZSBhdXRvY29ycmVsYXRpb24gaXMgc28gc3Ryb25nIGl0J3MgYXMgaWYgd2Ugb25seSBoYWQgNDY1IGluZGVwZW5kZW50IGRyYXdzIGZyb20gdGhlIHBvc3RlcmlvciBmb3Igc3RlZXBuZXNzISBUbyBpbXByb3ZlIHRoZSBlc3RpbWF0aW9uIG9mIHRoZXNlIHBhcmFtZXRlcnMsIG9uZSBjb3VsZCBydW4gdGhlIGNoYWlucyBsb25nZXIsIG9yIGNoYW5nZSB0aGUgbW9kZWwsIHN0cnVjdHVyYWxseSwgb3Igd2l0aCB0aGUgcHJpb3JzLgoKCgojIyBBcmUgdGhlcmUgY29ycmVsYXRpb25zIGJldHdlZW4gcGFyYW1ldGVycyBpbiBvdXIgbW9kZWw/CmBgYHtyfQpjcm9zc2NvcnIucGxvdChtY21jMikKYGBgCgpUaGVyZSBpcyBhIGNvcnJlbGF0aW9uIGJldHdlZW4gdHdvIHZhcmlhYmxlczoKYGBge3J9CmNyb3NzY29ycihtY21jMikKYGBgCgowLjY3OiBubyB3b25kZXIgc3RlZXBuZXNzIGFuZCBtaWRwb2ludCBhcmUgZGlmZmljdWx0IHRvIGVzdGltYXRlOiB0aGV5IGFyZSBzdHJvbmdseSBjb3JyZWxhdGVkLgoKIyMgQW5hbHlzaXMgb2YgdGhlIGVzdGltYXRlZCBwYXJhbWV0ZXJzCgpgYGB7cn0KCnBsb3QobWNtYzIpCmBgYAoKV2UgY2FuIHNlZSBvbiB0aGUgdHJhY2VzIHRoYXQgc3RlZXBuZXNzIGFuZCBtaWRwb2ludCBhcmUgZGlmZmljdWx0IHRvIHNhbXBsZSBmb3IgdGhlIG1jbWMuCgojIyBFc3RpbWF0ZWQgcGFyYW1ldGVycwpgYGB7cn0Kc3VtbWFyeShtY21jMikKYGBgCgoKCkxvb2tpbmcgYXQgdGhlIG1lZGlhbiB2YWx1ZXMsIHRoZSBlc3RpbWF0ZXMgYXJlIHByZXR0eSBnb29kIGZvciBtb3N0IHBhcmFtZXRlcnMuIFRoZSB0cnVlIHZhbHVlcyB3ZXJlIDEwMDAgZm9yIG1lYW5GcnVpdHMsIDUwMCBmb3IgbWlkcG9pbnQsIDUwIGZvciBvcHRpbWFsUHJlY2lwaXRhdGlvbiwgMC4xIGZvciBzdGVlcG5lc3MuIEJ1dCBpdCB3YXMgMTAgZm9yIHNjYWxlUHJpb3IsIHdoaWNoIGlzIGVzdGltYXRlZCBoZXJlIGF0IDEuMjcuIE9uZSB3b3VsZCBuZWVkIHRvIGludmVzdGlnYXRlIHdoeSBzY2FsZVByaW9yIGlzIG5vdCB3ZWxsIGVzdGltYXRlZC4KCg==