Commit fa2007cd authored by matbuoro's avatar matbuoro
Browse files

multiple chains included

parent 5bf98769
......@@ -44,7 +44,8 @@ if(!file.exists(paste("inits/init-",site,"-",stade,year,".txt",sep=""))){
#load(paste('inits/inits_',stade,'.Rdata',sep="")) # chargement des inits
#if(site == "Bresle" && stade == "adult") {inits <- list(read.bugsdata(paste("inits/init-",site,"-",stade,year,".txt",sep="")))}
#if(site == "Nivelle") {inits <- list(read.bugsdata(paste("inits/init-",site,"-",stade,year,".txt",sep="")))}
inits <- list(read.bugsdata(paste("inits/init-",site,"-",stade,year,".txt",sep="")))
inits.tmp <- read.bugsdata(paste("inits/init-",site,"-",stade,year,".txt",sep=""))
inits <- rep(list(inits.tmp),3)
#------------------------MODEL----------------------------------##
model <- paste("model/model_",stade,"-",site,".R",sep="") # path of the model
......@@ -56,11 +57,11 @@ filename <- file.path(work.dir, model)
#---------------------------ANALYSIS-----------------------------##
nChains = length(inits) # Number of chains to run.
nChains = 3 #length(inits) # Number of chains to run.
adaptSteps = 1000 # Number of steps to "tune" the samplers.
nburnin=5000 # Number of steps to "burn-in" the samplers.
nstore=25000 # Total number of steps in chains to save.
nthin=2 # Number of steps to "thin" (1=keep every step).
nthin=4 # Number of steps to "thin" (1=keep every step).
#nPerChain = ceiling( ( numSavedSteps * thinSteps ) / nChains ) # Steps per chain.
### Start of the run ###
......@@ -139,7 +140,7 @@ cat("=============================\n")
if (nChains > 1) {
cat("Convergence: gelman-Rubin R test\n")
gelman.diag(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)])
gelman.diag(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)],multivariate=TRUE)
}
......@@ -177,6 +178,7 @@ pdf(paste('results/Results_',stade,"_",year,'.pdf',sep=""))
traplot(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)])
# caterplot(fit.mcmc,parameters[i])
#}
gelman.plot(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)])
dev.off()
......
OpenBUGS did not run correctly.
OpenBUGS did not run correctly.
OpenBUGS did not run correctly.
list(lambda_tot0=1.04200E+02, logit_flow_Eu=c(-4.26200E-01, -7.07900E-01), lflow_fall_Eu=c(-5.00000E-01, -5.00000E-01), logit_int_Eu=c(5.08200E-01, 8.30400E-01), mupi_B=c(8.27000E-02, 1.03200E-01), pi_Eu00=c(5.87100E-01, 5.99200E-01), pi_Eu01=c(2.00600E-01, 7.99400E-01), rate_lambda=4.23800E-02, s=c(1.56700E+01, 3.91400E+00), shape_lambda=5.99500E+00, sigmapi_B=c(6.44000E-01, 6.96900E-01), sigmapi_Eu=c(1.01900E+00, 8.86500E-01), lambda_tot=c(1.02000E+02, 1.15500E+02, 2.20800E+02, 1.76200E+02, 1.30000E+02, NA, 1.24500E+02, 1.94300E+02, 1.57000E+02, NA, 9.30000E+01, 1.44000E+02, 6.23300E+01, 3.40000E+01, 2.80500E+02, 1.75000E+02, NA, NA, 1.30300E+02, 2.30000E+01, 8.72300E+01, 3.18200E+02, 2.43200E+02, 2.64000E+02, 1.62000E+02, 1.09400E+02, 1.93500E+02, 1.47000E+02, 1.51200E+02, 2.97000E+02, 2.02400E+02, 1.58900E+02, 1.96500E+02), logit_pi_Eu= structure(.Data= c(1.68200E+00, 3.92200E-02, NA, -7.12400E-01, NA, -1.57100E+00, NA, -9.26000E-01, 1.48500E+00, -1.48500E+00, NA, NA, 1.40600E+00, 1.20600E-01, NA, -1.41400E+00, NA, -2.17400E-01, NA, NA, 2.37700E-01, -2.11600E+00, -7.56300E-01, -4.26300E+00, 3.12700E-01, 5.34900E-02, NA, 3.56700E-01, NA, -2.56700E+00, -1.83900E+00, -1.74600E+00, NA, NA, NA, NA, 1.81000E+00, -1.38900E+00, NA, 1.89700E+00, NA, -1.34700E+00, 7.63400E+00, -1.45900E+00, 4.06300E-01, 1.71000E-01, -1.18100E+00, -2.21400E+00, 3.25800E+00, -3.25800E+00, NA, 4.95300E-01, NA, -1.10900E+00, NA, -7.23900E-01, NA, -7.66000E-01, 5.16200E-01, -9.30500E-01, 5.42600E-01, -9.61000E-01, NA, -3.41400E-01, NA, -1.71400E+00), .Dim=c(33, 2)), logit_pi_B= structure(.Data= c(-2.46400E+00, -4.61500E+00, -1.33800E+00, -4.74000E+00, -2.69300E+00, -4.69500E+00, -2.92200E+00, -3.76200E+00, -4.86000E+00, -4.86000E+00, NA, NA, -3.70100E+00, NA, -2.91400E+00, NA, -3.41400E+00, -4.35000E+00, NA, NA, -3.81800E+00, -4.52200E+00, -2.01100E+00, -4.96300E+00, -2.24000E+00, -2.43900E+00, -3.49700E+00, NA, -3.10800E+00, -4.52700E+00, -4.46000E+00, -3.52600E+00, NA, NA, NA, NA, -1.43800E+00, -2.60100E+00, NA, -3.09100E+00, -1.34700E+00, -3.33500E+00, -1.43900E+00, -4.13700E+00, -2.46800E+00, -3.25900E+00, -3.76100E+00, NA, -2.83300E+00, -5.08100E+00, -2.00400E+00, -3.27200E+00, -2.63100E+00, -3.63000E+00, -2.17500E+00, -4.98400E+00, -2.76000E+00, -5.01200E+00, -3.58700E+00, -4.58500E+00, -3.48800E+00, -3.67600E+00, -2.50500E+00, -3.65600E+00, -2.82500E+00, -5.27600E+00), .Dim=c(33, 2)), n= structure(.Data= c(6.40000E+01, 3.90000E+01, 9.30000E+01, 2.30000E+01, 1.91000E+02, 3.10000E+01, 1.41000E+02, 3.60000E+01, 1.06000E+02, 2.40000E+01, NA, NA, 7.50000E+01, 5.00000E+01, 1.68000E+02, 2.70000E+01, 1.23000E+02, 3.50000E+01, NA, NA, 7.80000E+01, 1.50000E+01, 1.38000E+02, 6.00000E+00, 3.30000E+01, 3.00000E+01, 2.40000E+01, 1.00000E+01, 2.66000E+02, 1.50000E+01, 8.40000E+01, 9.10000E+01, NA, NA, NA, NA, 1.06000E+02, 2.50000E+01, 1.30000E+01, 1.00000E+01, 7.30000E+01, 1.50000E+01, 2.68000E+02, 5.10000E+01, 1.28000E+02, 1.16000E+02, 1.86000E+02, 7.80000E+01, 1.56000E+02, 6.00000E+00, 7.50000E+01, 3.60000E+01, 1.58000E+02, 3.60000E+01, 1.23000E+02, 2.40000E+01, 1.22000E+02, 3.00000E+01, 2.05000E+02, 9.30000E+01, 1.41000E+02, 6.20000E+01, 1.22000E+02, 3.80000E+01, 1.74000E+02, 2.30000E+01), .Dim=c(33, 2)))
list(lambda_tot0=1.04200E+02, logit_flow_Eu=c(-4.26200E-01, -7.07900E-01), lflow_fall_Eu=c(-5.00000E-01, -5.00000E-01), logit_int_Eu=c(5.08200E-01, 8.30400E-01), mupi_B=c(8.27000E-02, 1.03200E-01), pi_Eu00=c(5.87100E-01, 5.99200E-01), pi_Eu01=c(2.00600E-01, 7.99400E-01), rate_lambda=4.23800E-02, s=c(1.56700E+01, 3.91400E+00), shape_lambda=5.99500E+00, sigmapi_B=c(6.44000E-01, 6.96900E-01), sigmapi_Eu=c(1.01900E+00, 8.86500E-01), lambda_tot=c(1.02000E+02, 1.15500E+02, 2.20800E+02, 1.76200E+02, 1.30000E+02, NA, 1.24500E+02, 1.94300E+02, 1.57000E+02, NA, 9.30000E+01, 1.44000E+02, 6.23300E+01, 3.40000E+01, 2.80500E+02, 1.75000E+02, NA, NA, 1.30300E+02, 2.30000E+01, 8.72300E+01, 3.18200E+02, 2.43200E+02, 2.64000E+02, 1.62000E+02, 1.09400E+02, 1.93500E+02, 1.47000E+02, 1.51200E+02, 2.97000E+02, 2.02400E+02, 1.58900E+02, 1.96500E+02), logit_pi_Eu= structure(.Data= c(1.68200E+00, 3.92200E-02, NA, -7.12400E-01, NA, -1.57100E+00, NA, -9.26000E-01, 1.48500E+00, -1.48500E+00, NA, NA, 1.40600E+00, 1.20600E-01, NA, -1.41400E+00, NA, -2.17400E-01, NA, NA, 2.37700E-01, -2.11600E+00, -7.56300E-01, -4.26300E+00, 3.12700E-01, 5.34900E-02, NA, 3.56700E-01, NA, -2.56700E+00, -1.83900E+00, -1.74600E+00, NA, NA, NA, NA, 1.81000E+00, -1.38900E+00, NA, 1.89700E+00, NA, -1.34700E+00, 7.63400E+00, -1.45900E+00, 4.06300E-01, 1.71000E-01, -1.18100E+00, -2.21400E+00, 3.25800E+00, -3.25800E+00, NA, 4.95300E-01, NA, -1.10900E+00, NA, -7.23900E-01, NA, -7.66000E-01, 5.16200E-01, -9.30500E-01, 5.42600E-01, -9.61000E-01, NA, -3.41400E-01, NA, -1.71400E+00), .Dim=c(33, 2)), logit_pi_B= structure(.Data= c(-2.46400E+00, -4.61500E+00, -1.33800E+00, -4.74000E+00, -2.69300E+00, -4.69500E+00, -2.92200E+00, -3.76200E+00, -4.86000E+00, -4.86000E+00, NA, NA, -3.70100E+00, NA, -2.91400E+00, NA, -3.41400E+00, -4.35000E+00, NA, NA, -3.81800E+00, -4.52200E+00, -2.01100E+00, -4.96300E+00, -2.24000E+00, -2.43900E+00, -3.49700E+00, NA, -3.10800E+00, -4.52700E+00, -4.46000E+00, -3.52600E+00, NA, NA, NA, NA, -1.43800E+00, -2.60100E+00, NA, -3.09100E+00, -1.34700E+00, -3.33500E+00, -1.43900E+00, -4.13700E+00, -2.46800E+00, -3.25900E+00, -3.76100E+00, NA, -2.83300E+00, -5.08100E+00, -2.00400E+00, -3.27200E+00, -2.63100E+00, -3.63000E+00, -2.17500E+00, -4.98400E+00, -2.76000E+00, -5.01200E+00, -3.58700E+00, -4.58500E+00, -3.48800E+00, -3.67600E+00, -2.50500E+00, -3.65600E+00, -2.82500E+00, -5.27600E+00), .Dim=c(33, 2)), n= structure(.Data= c(6.40000E+01, 3.90000E+01, 9.30000E+01, 2.30000E+01, 1.91000E+02, 3.10000E+01, 1.41000E+02, 3.60000E+01, 1.06000E+02, 2.40000E+01, NA, NA, 7.50000E+01, 5.00000E+01, 1.68000E+02, 2.70000E+01, 1.23000E+02, 3.50000E+01, NA, NA, 7.80000E+01, 1.50000E+01, 1.38000E+02, 6.00000E+00, 3.30000E+01, 3.00000E+01, 2.40000E+01, 1.00000E+01, 2.66000E+02, 1.50000E+01, 8.40000E+01, 9.10000E+01, NA, NA, NA, NA, 1.06000E+02, 2.50000E+01, 1.30000E+01, 1.00000E+01, 7.30000E+01, 1.50000E+01, 2.68000E+02, 5.10000E+01, 1.28000E+02, 1.16000E+02, 1.86000E+02, 7.80000E+01, 1.56000E+02, 6.00000E+00, 7.50000E+01, 3.60000E+01, 1.58000E+02, 3.60000E+01, 1.23000E+02, 2.40000E+01, 1.22000E+02, 3.00000E+01, 2.05000E+02, 9.30000E+01, 1.41000E+02, 6.20000E+01, 1.22000E+02, 3.80000E+01, 1.74000E+02, 2.30000E+01), .Dim=c(33, 2)))
This diff is collapsed.
modelCheck('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Bresle/adult/bugs/model_adult-Bresle.R.txt')
modelData('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Bresle/adult/bugs/data.txt')
modelCompile(1)
modelCompile(3)
modelSetRN(1)
modelInits('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Bresle/adult/bugs/inits1.txt',1)
modelInits('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Bresle/adult/bugs/inits2.txt',2)
modelInits('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Bresle/adult/bugs/inits3.txt',3)
modelGenInits()
modelUpdate(5000,2,5000)
modelUpdate(5000,4,5000)
samplesSet(pi_Eu00)
samplesSet(pi_Eu01)
samplesSet(logit_int_Eu)
......@@ -55,7 +57,7 @@ summarySet(Plambda0)
summarySet(s)
summarySet(lambda_tot)
summarySet(Plambda)
modelUpdate(25000,2,25000)
modelUpdate(25000,4,25000)
samplesCoda('*', '/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Bresle/adult/bugs//')
summaryStats('*')
modelQuit('y')
......@@ -12,15 +12,15 @@ heidel.diag also implements a convergence diagnostic, and removes up to half the
Stationarity start p-value
test iteration
shape_lambda passed 10001 0.223
rate_lambda passed 2501 0.157
lambda_tot0 passed 1 0.593
Halfwidth Mean Halfwidth
test
shape_lambda passed 3.627 0.06632
rate_lambda passed 0.024 0.00042
lambda_tot0 passed 165.669 1.74561
shape_lambda passed 1 0.7430
rate_lambda passed 1 0.3791
lambda_tot0 passed 5001 0.0613
Halfwidth Mean Halfwidth
test
shape_lambda passed 3.6844 0.047101
rate_lambda passed 0.0243 0.000303
lambda_tot0 passed 168.3312 1.382674
---------------------------
Geweke's convergence diagnostic
......@@ -36,7 +36,7 @@ Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5
shape_lambda rate_lambda lambda_tot0
2.851 3.231 -0.877
0.910 0.476 -1.454
---------------------------
......@@ -48,7 +48,7 @@ Probability (s) = 0.95
Burn-in Total Lower bound Dependence
(M) (N) (Nmin) factor (I)
shape_lambda 24 30996 3746 8.27
rate_lambda 24 30240 3746 8.07
lambda_tot0 12 17600 3746 4.70
shape_lambda 24 34744 3746 9.27
rate_lambda 32 35472 3746 9.47
lambda_tot0 12 18352 3746 4.90
This diff is collapsed.
......@@ -44,7 +44,8 @@ if(!file.exists(paste("inits/init-",site,"-",stade,year,".txt",sep=""))){
#load(paste('inits/inits_',stade,'.Rdata',sep="")) # chargement des inits
#if(site == "Bresle" && stade == "adult") {inits <- list(read.bugsdata(paste("inits/init-",site,"-",stade,year,".txt",sep="")))}
#if(site == "Nivelle") {inits <- list(read.bugsdata(paste("inits/init-",site,"-",stade,year,".txt",sep="")))}
inits <- list(read.bugsdata(paste("inits/init-",site,"-",stade,year,".txt",sep="")))
inits.tmp <- read.bugsdata(paste("inits/init-",site,"-",stade,year,".txt",sep=""))
inits <- rep(list(inits.tmp),3)
#------------------------MODEL----------------------------------##
model <- paste("model/model_",stade,"-",site,".R",sep="") # path of the model
......@@ -56,11 +57,11 @@ filename <- file.path(work.dir, model)
#---------------------------ANALYSIS-----------------------------##
nChains = length(inits) # Number of chains to run.
nChains = 3 #length(inits) # Number of chains to run.
adaptSteps = 1000 # Number of steps to "tune" the samplers.
nburnin=5000 # Number of steps to "burn-in" the samplers.
nstore=25000 # Total number of steps in chains to save.
nthin=2 # Number of steps to "thin" (1=keep every step).
nthin=4 # Number of steps to "thin" (1=keep every step).
#nPerChain = ceiling( ( numSavedSteps * thinSteps ) / nChains ) # Steps per chain.
### Start of the run ###
......@@ -139,7 +140,7 @@ cat("=============================\n")
if (nChains > 1) {
cat("Convergence: gelman-Rubin R test\n")
gelman.diag(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)])
gelman.diag(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)],multivariate=TRUE)
}
......@@ -177,6 +178,7 @@ pdf(paste('results/Results_',stade,"_",year,'.pdf',sep=""))
traplot(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)])
# caterplot(fit.mcmc,parameters[i])
#}
gelman.plot(fit.mcmc[,which(varnames(fit.mcmc)%in%parameterstotest)])
dev.off()
......
list(mu_B=5.00000E-01, sigmap_B=1.00000E+00, logit_int_Eu=1.00000E+00, logit_flow_Eu=1.00000E+00, sigmap_Eu=1.00000E+00, shape_lambda=2.50000E+00, rate_lambda=1.00000E-02, Ntot=c( NA, NA, NA, NA, 1.93100E+03, NA, NA, NA, NA, NA, NA, 1.65000E+03, 2.36600E+03, NA, 3.10000E+03, 6.26700E+03, 1.59900E+03, 6.66000E+02, 1.63200E+03, NA, NA, 2.91700E+03, 7.56700E+03, 5.21500E+03, 2.72000E+03, 5.54400E+03, 6.58700E+03, 1.05500E+03, 3.09100E+03, 6.40600E+03, 6.37300E+03, 1.97300E+03, 3.02000E+03, 7.85000E+03, 7.84200E+03), lambda=c( NA, NA, NA, NA, 1.92100E+03, NA, NA, NA, NA, NA, NA, 1.64000E+03, 2.35600E+03, NA, 3.09000E+03, 6.25700E+03, 1.58900E+03, 6.56000E+02, 1.62200E+03, NA, NA, 2.90700E+03, 7.55700E+03, 5.20500E+03, 2.71000E+03, 5.53400E+03, 6.57700E+03, 1.04500E+03, 3.08100E+03, 6.39600E+03, 6.36300E+03, 1.96300E+03, 3.01000E+03, 7.84000E+03, 7.83200E+03), logit_pi_B=c(5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01), logit_pi_Eu=c(5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01))
list(mu_B=5.00000E-01, sigmap_B=1.00000E+00, logit_int_Eu=1.00000E+00, logit_flow_Eu=1.00000E+00, sigmap_Eu=1.00000E+00, shape_lambda=2.50000E+00, rate_lambda=1.00000E-02, Ntot=c( NA, NA, NA, NA, 1.93100E+03, NA, NA, NA, NA, NA, NA, 1.65000E+03, 2.36600E+03, NA, 3.10000E+03, 6.26700E+03, 1.59900E+03, 6.66000E+02, 1.63200E+03, NA, NA, 2.91700E+03, 7.56700E+03, 5.21500E+03, 2.72000E+03, 5.54400E+03, 6.58700E+03, 1.05500E+03, 3.09100E+03, 6.40600E+03, 6.37300E+03, 1.97300E+03, 3.02000E+03, 7.85000E+03, 7.84200E+03), lambda=c( NA, NA, NA, NA, 1.92100E+03, NA, NA, NA, NA, NA, NA, 1.64000E+03, 2.35600E+03, NA, 3.09000E+03, 6.25700E+03, 1.58900E+03, 6.56000E+02, 1.62200E+03, NA, NA, 2.90700E+03, 7.55700E+03, 5.20500E+03, 2.71000E+03, 5.53400E+03, 6.57700E+03, 1.04500E+03, 3.08100E+03, 6.39600E+03, 6.36300E+03, 1.96300E+03, 3.01000E+03, 7.84000E+03, 7.83200E+03), logit_pi_B=c(5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01), logit_pi_Eu=c(5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01, 5.00000E-01))
This diff is collapsed.
modelCheck('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Bresle/smolt/bugs/model_smolt-Bresle.R.txt')
modelData('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Bresle/smolt/bugs/data.txt')
modelCompile(1)
modelCompile(3)
modelSetRN(1)
modelInits('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Bresle/smolt/bugs/inits1.txt',1)
modelInits('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Bresle/smolt/bugs/inits2.txt',2)
modelInits('/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Bresle/smolt/bugs/inits3.txt',3)
modelGenInits()
modelUpdate(5000,2,5000)
modelUpdate(5000,4,5000)
samplesSet(mu_B)
samplesSet(sigmap_B)
samplesSet(logit_int_Eu)
......@@ -49,7 +51,7 @@ summarySet(lambda)
summarySet(Nesc)
summarySet(test)
summarySet(R2)
modelUpdate(25000,2,25000)
modelUpdate(25000,4,25000)
samplesCoda('*', '/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance/Bresle/smolt/bugs//')
summaryStats('*')
modelQuit('y')
=============================
DIAGNOSTICS
=============================
Convergence: gelman-Rubin R test
Potential scale reduction factors:
Point est. Upper C.I.
mu_B 1.00 1.01
sigmap_B 1.01 1.02
logit_int_Eu 1.00 1.00
logit_flow_Eu 1.00 1.01
sigmap_Eu 1.00 1.00
p_B95 1.00 1.01
p_B00 1.02 1.07
p_B02 1.00 1.01
shape_lambda 1.01 1.02
rate_lambda 1.01 1.02
mean_gamma 1.00 1.01
var_gamma 1.01 1.02
test 1.01 1.01
R2 1.01 1.02
Multivariate psrf
1.03
---------------------------
Heidelberger and Welch's convergence diagnostic
......@@ -9,40 +31,114 @@ heidel.diag is a run length control diagnostic based on a criterion of relative
heidel.diag also implements a convergence diagnostic, and removes up to half the chain in order to ensure that the means are estimated from a chain that has converged.
[[1]]
Stationarity start p-value
test iteration
mu_B passed 10001 0.087859
sigmap_B passed 1 0.055910
logit_int_Eu passed 1 0.134761
logit_flow_Eu passed 10001 0.268373
sigmap_Eu passed 7501 0.080705
p_B95 passed 1 0.356449
p_B00 passed 1 0.160608
p_B02 passed 1 0.444178
shape_lambda passed 5001 0.084974
rate_lambda passed 1 0.115752
mean_gamma failed NA 0.000959
var_gamma passed 10001 0.187606
test passed 10001 0.116361
R2 passed 10001 0.083408
mu_B failed NA 5.62e-07
sigmap_B passed 1 1.50e-01
logit_int_Eu passed 1 4.55e-01
logit_flow_Eu passed 1 5.75e-01
sigmap_Eu passed 1 9.25e-01
p_B95 passed 1 6.74e-01
p_B00 failed NA 8.71e-03
p_B02 passed 1 1.88e-01
shape_lambda failed NA 3.38e-04
rate_lambda failed NA 1.66e-02
mean_gamma failed NA 2.15e-06
var_gamma passed 7501 1.26e-01
test passed 1 4.72e-01
R2 passed 1 8.25e-01
Halfwidth Mean Halfwidth
test
mu_B passed 2.74e-01 3.50e-03
sigmap_B passed 8.71e-01 1.49e-02
logit_int_Eu passed -2.50e+00 4.00e-03
logit_flow_Eu passed -7.42e-02 5.34e-03
sigmap_Eu passed 2.95e-01 4.84e-03
p_B95 failed 2.45e-01 4.53e-02
p_B00 passed 6.12e-01 3.73e-02
p_B02 failed 2.30e-01 4.48e-02
shape_lambda passed 2.61e+00 5.39e-02
rate_lambda passed 6.85e-04 1.09e-05
mu_B <NA> NA NA
sigmap_B passed 9.05e-01 1.93e-02
logit_int_Eu passed -2.50e+00 2.84e-03
logit_flow_Eu passed -7.47e-02 3.85e-03
sigmap_Eu passed 2.95e-01 3.04e-03
p_B95 failed 2.28e-01 3.89e-02
p_B00 <NA> NA NA
p_B02 failed 2.86e-01 2.98e-02
shape_lambda <NA> NA NA
rate_lambda <NA> NA NA
mean_gamma <NA> NA NA
var_gamma passed 6.32e+06 1.28e+05
test passed 1.54e-01 1.51e-02
R2 failed 2.61e-02 8.95e-03
var_gamma passed 6.04e+06 1.17e+05
test passed 1.58e-01 1.36e-02
R2 failed 2.83e-02 6.54e-03
[[2]]
Stationarity start p-value
test iteration
mu_B passed 10001 5.35e-02
sigmap_B failed NA 1.05e-06
logit_int_Eu passed 10001 4.40e-01
logit_flow_Eu passed 1 3.63e-01
sigmap_Eu passed 1 2.56e-01
p_B95 passed 1 4.46e-01
p_B00 passed 1 2.86e-01
p_B02 passed 1 6.67e-01
shape_lambda passed 1 6.46e-02
rate_lambda passed 5001 8.84e-02
mean_gamma passed 2501 8.42e-02
var_gamma passed 2501 5.30e-02
test passed 1 3.80e-01
R2 passed 1 6.06e-01
Halfwidth Mean Halfwidth
test
mu_B passed 2.94e-01 3.63e-03
sigmap_B <NA> NA NA
logit_int_Eu passed -2.50e+00 3.84e-03
logit_flow_Eu passed -7.13e-02 3.85e-03
sigmap_Eu passed 2.96e-01 2.94e-03
p_B95 failed 2.51e-01 3.56e-02
p_B00 passed 5.63e-01 4.49e-02
p_B02 failed 2.74e-01 3.70e-02
shape_lambda passed 2.48e+00 5.67e-02
rate_lambda passed 6.70e-04 1.72e-05
mean_gamma passed 3.74e+03 3.07e+01
var_gamma passed 6.22e+06 1.71e+05
test passed 1.67e-01 1.46e-02
R2 failed 2.58e-02 6.17e-03
[[3]]
Stationarity start p-value
test iteration
mu_B failed NA 0.000197
sigmap_B passed 1 0.543246
logit_int_Eu passed 1 0.107623
logit_flow_Eu passed 1 0.122057
sigmap_Eu passed 1 0.113756
p_B95 passed 1 0.588786
p_B00 passed 1 0.278269
p_B02 passed 1 0.129660
shape_lambda failed NA 0.001295
rate_lambda passed 1 0.079150
mean_gamma passed 2501 0.060192
var_gamma passed 1 0.157665
test passed 10001 0.269950
R2 passed 1 0.098299
Halfwidth Mean Halfwidth
test
mu_B <NA> NA NA
sigmap_B passed 9.07e-01 1.60e-02
logit_int_Eu passed -2.50e+00 2.76e-03
logit_flow_Eu passed -7.94e-02 3.60e-03
sigmap_Eu passed 2.95e-01 3.05e-03
p_B95 failed 2.20e-01 3.76e-02
p_B00 passed 5.87e-01 4.20e-02
p_B02 failed 2.73e-01 3.51e-02
shape_lambda <NA> NA NA
rate_lambda passed 6.60e-04 1.49e-05
mean_gamma passed 3.77e+03 3.82e+01
var_gamma passed 6.39e+06 1.54e+05
test passed 1.28e-01 1.11e-02
R2 failed 3.88e-02 5.82e-03
---------------------------
Geweke's convergence diagnostic
......@@ -53,20 +149,97 @@ The test statistic is a standard Z-score: the difference between the two sample
The Z-score is calculated under the assumption that the two parts of the chain are asymptotically independent, which requires that the sum of frac1 and frac2 be strictly less than 1.
[[1]]
Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5
mu_B sigmap_B logit_int_Eu logit_flow_Eu sigmap_Eu
0.436 -1.226 0.253 0.553 0.371
p_B95 p_B00 p_B02 shape_lambda rate_lambda
0.248 -0.610 0.897 -1.056 -0.701
mean_gamma var_gamma test R2
-2.617 -0.339 0.495 -0.867
[[2]]
Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5
mu_B sigmap_B logit_int_Eu logit_flow_Eu sigmap_Eu
2.439 2.125 -1.220 -2.111 0.848
p_B95 p_B00 p_B02 shape_lambda rate_lambda
1.029 0.582 0.516 0.504 2.145
mean_gamma var_gamma test R2
-4.928 -3.915 -2.464 1.559
[[3]]
Fraction in 1st window = 0.1
Fraction in 2nd window = 0.5
mu_B sigmap_B logit_int_Eu logit_flow_Eu sigmap_Eu
7.863 -3.134 2.159 -2.086 -1.054
0.612 -1.158 1.042 -0.236 0.909
p_B95 p_B00 p_B02 shape_lambda rate_lambda
-1.072 -2.978 -1.147 -2.020 0.354
-0.499 -1.550 1.834 0.323 1.307
mean_gamma var_gamma test R2
-10.903 -2.704 -2.770 3.015
-1.997 -2.503 0.396 0.322
---------------------------
Raftery and Lewis's diagnostic
[[1]]
Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95
Burn-in Total Lower bound Dependence
(M) (N) (Nmin) factor (I)
mu_B 48 78128 3746 20.90
sigmap_B 36 51216 3746 13.70
logit_int_Eu 8 16040 3746 4.28
logit_flow_Eu 36 47808 3746 12.80
sigmap_Eu 36 48144 3746 12.90
p_B95 320 319008 3746 85.20
p_B00 260 328224 3746 87.60
p_B02 512 617984 3746 165.00
shape_lambda 24 35120 3746 9.38
rate_lambda 16 19520 3746 5.21
mean_gamma 8 15840 3746 4.23
var_gamma 24 36272 3746 9.68
test 96 856448 3746 229.00
R2 8 15180 3746 4.05
[[2]]
Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
Probability (s) = 0.95
Burn-in Total Lower bound Dependence
(M) (N) (Nmin) factor (I)
mu_B 12 16180 3746 4.32
sigmap_B 48 66128 3746 17.70
logit_int_Eu 12 16608 3746 4.43
logit_flow_Eu 24 31064 3746 8.29
sigmap_Eu 60 85240 3746 22.80
p_B95 884 902304 3746 241.00
p_B00 384 426960 3746 114.00
p_B02 312 396500 3746 106.00
shape_lambda 60 84620 3746 22.60
rate_lambda 16 20840 3746 5.56
mean_gamma 8 15596 3746 4.16
var_gamma 24 32872 3746 8.78
test 128 859328 3746 229.00
R2 8 15308 3746 4.09
[[3]]
Quantile (q) = 0.025
Accuracy (r) = +/- 0.005
......@@ -74,18 +247,19 @@ Probability (s) = 0.95
Burn-in Total Lower bound Dependence
(M) (N) (Nmin) factor (I)
mu_B 12 21996 3746 5.87
sigmap_B 18 25692 3746 6.86
logit_int_Eu 12 17680 3746 4.72
logit_flow_Eu 8 15736 3746 4.20
sigmap_Eu 18 25518 3746 6.81
p_B95 266 294120 3746 78.50
p_B00 264 303908 3746 81.10
p_B02 468 487578 3746 130.00
shape_lambda 20 20072 3746 5.36
rate_lambda 16 17528 3746 4.68
mean_gamma 4 7766 3746 2.07
var_gamma 24 28434 3746 7.59
test 30 220960 3746 59.00
R2 4 7832 3746 2.09
mu_B 12 16468 3746 4.40
sigmap_B 64 69280 3746 18.50
logit_int_Eu 24 31488 3746 8.41
logit_flow_Eu 24 47088 3746 12.60
sigmap_Eu 80 91820 3746 24.50
p_B95 416 520780 3746 139.00
p_B00 256 271680 3746 72.50
p_B02 572 662948 3746 177.00
shape_lambda 24 37264 3746 9.95
rate_lambda 24 33232 3746 8.87
mean_gamma 8 15832 3746 4.23
var_gamma 24 36152 3746 9.65
test 72 510888 3746 136.00
R2 8 14816 3746 3.96
This diff is collapsed.
......@@ -6,7 +6,7 @@ YEAR=2016
CHAINS=3
BURNIN=5000 # Number of steps to "burn-in" the samplers.
ITER=25000 # Total number of steps in chains to save.
THIN=2 # Number of steps to "thin" (1=keep every step).
THIN=4 # Number of steps to "thin" (1=keep every step).
REPbase="/home/basp-meco88/Documents/RESEARCH/PROJECTS/ORE/Abundance"
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment