accrate.depth()
does only return non NA
values for MCMC samples when an input age is not reached. This means that the output vector length will change, depending on how often an age value is not within the MCMC samples. In addition, the order of the returned values depends on the depth intervals the defined age is within.
Intuitively, I would expect that the order of the MCMC samples should not be changed since these are conditionally not independent and I would expect the vector to alwas have the full length of potential MCMC samples with NA
s where ages are not within the MCMC sample range.
Below, there is a reproducible example for this behavior and a modified function accrate.depth_new()
which returns accumulation rates in a format I would expect (the same format also other functions, e.g. Bacon.d.Age()
or Bacon.Age.d()
return).
library(rbacon)
#> Warning: package 'rbacon' was built under R version 4.0.5
#> Loading required package: IntCal
#> Warning: package 'IntCal' was built under R version 4.0.5
#compute age-depth model for demo core
rbacon::Bacon(coredir = getwd(), ask = FALSE, accept.suggestions = TRUE, verbose = FALSE, plot.pdf = FALSE, seed = 2312)
#> The run's files will be put in this folder: x\AppData\Local\Temp\RtmpaiJGtg\reprex-4dac6510213e-fast-perch/MSB2K
#> Will run 4646000 iterations and store around 2019
#> Reading x\AppData\Local\Temp\RtmpaiJGtg\reprex-4dac6510213e-fast-perch/MSB2K/MSB2K_21.bacon
#> Constant calibration curve.
#> IntCal20: Reading from file: x\Documents\R\win-library\4.0\IntCal\extdata/3Col_intcal20.14C
#> Marine20: Reading from file: x\R\win-library\4.0\IntCal\extdata/3Col_marine20.14C
#> SHCal20: Reading from file: x\Documents\R\win-library\4.0\IntCal\extdata/3Col_shcal20.14C
#> Added det: GrA-19478: 4128.0+-65.0 d=1.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19143: 4106.0+-60.0 d=4.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19144: 4046.0+-59.0 d=8.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19146: 4184.0+-58.0 d=12.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19147: 4076.0+-62.0 d=14.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19148: 4107.0+-61.0 d=14.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19141: 4097.0+-58.0 d=14.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-18675: 4177.0+-53.0 d=17.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19151: 4220.0+-59.0 d=20.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19476: 4281.0+-64.0 d=21.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19475: 4374.0+-64.0 d=21.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19152: 4493.0+-62.0 d=22.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-18323: 4452.0+-52.0 d=28.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19474: 4616.0+-64.0 d=31.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19473: 4662.0+-64.0 d=32.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19509: 4743.0+-67.0 d=33.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19480: 4638.0+-67.0 d=34.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19483: 4810.0+-67.0 d=37.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-18674: 4757.0+-82.0 d=38.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19153: 4839.0+-59.0 d=41.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19484: 4913.0+-65.0 d=43.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19154: 4880.0+-57.0 d=46.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19485: 4989.0+-70.0 d=47.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19486: 5070.0+-66.0 d=48.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19488: 4993.0+-67.0 d=49.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-17626: 5115.0+-79.0 d=50.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-18682: 5026.0+-51.0 d=52.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19489: 5242.0+-64.0 d=53.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-18679: 5159.0+-50.0 d=54.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19490: 5130.0+-66.0 d=55.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19492: 5238.0+-65.0 d=58.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-17501: 5293.0+-38.0 d=59.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-18320: 5293.0+-54.0 d=64.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-18678: 5368.0+-51.0 d=70.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-19494: 5498.0+-69.0 d=71.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-18319: 5588.0+-55.0 d=73.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-18318: 5514.0+-57.0 d=75.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-18688: 5535.0+-52.0 d=77.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-17627: 5644.0+-77.0 d=79.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> Added det: GrA-17508: 5885.0+-45.0 d=99.5 ResCorr= 0.0+-0.0 a=3 b=4 cc=IntCal20
#> BaconFixed: Bacon jumps model with fixed c's.
#> K=21, H=0, dim=23, Seed=2312, Dc=5.000000, c(0)=1.500000, c(K)=106.500000
#>
#> twalk: 5060000 iterations to run, (nil)
#> twalk thinning: 1 out of every 115 accepted iterations will be saved in file x\AppData\Local\Temp\RtmpaiJGtg\reprex-4dac6510213e-fast-perch/MSB2K/MSB2K_21.out
#> twalk: Finished, 0.7% of moved pars per iteration (ratio 35679.217391/5060000). Output in file (nil),
#> x\AppData\Local\Temp\RtmpaiJGtg\reprex-4dac6510213e-fast-perch/MSB2K/MSB2K_21.out
#> bacon: burn in (initial iterations which will be removed): 23000
#> Eso es to...eso es to...eso es to...eso es toooodo amigos!
#> Calculating age ranges...
#> | | | 0% | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |========= | 12% | |========= | 13% | |========== | 14% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============== | 19% | |============== | 20% | |=============== | 21% | |================ | 22% | |================ | 23% | |================= | 24% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================= | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |============================================ | 62% | |============================================ | 63% | |============================================= | 64% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================= | 69% | |================================================= | 70% | |================================================== | 71% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%
#>
#> Preparing ghost graph...
#> | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%
#>
#> Mean 95% confidence ranges 207 yr, min. 111 yr at 73.5 cm, max. 321 yr at 25.5 cm
#> 100% of the dates overlap with the age-depth model (95% ranges)
#>
dev.off()
#> null device
#> 1
x11()
# modified function accrate.depth()
accrate.age_new <- function(age, set=get('info'), cmyr=FALSE, ages=c(), BCAD=set$BCAD, silent=TRUE) {
if(length(ages) == 0) {
ages <- array(0, dim=c(nrow(set$output), length(set$elbows)))
for(i in 1:ncol(ages))
ages[,i] <- Bacon.Age.d(set$elbows[i], BCAD=BCAD)
}
if(!silent)
if(age < min(ages) || age > max(ages))
stop(" Warning, age outside the core's age range!\n")
# can this be made faster, i.e. without the loop?
accs <- rep(NA_real_, nrow(ages))
for(i in 2:ncol(ages)) {
if(BCAD)
these <- (ages[,i-1] > age) * (ages[,i] < age) else
these <- (ages[,i-1] < age) * (ages[,i] > age)
if(sum(these) > 0) # age lies within these age-model iterations
accs[these>0] <- set$output[which(these>0),i] # was i+1
}
if(cmyr)
accs <- 1/accs
return(accs)
}
# compute accumulation rate for an age where there are no NAs in MCMC samples
a1 <- rbacon::accrate.age(5000)
a2 <- accrate.age_new(5000)
plot(a1 ~ a2) # would expect all points on the identity line
hist(a1) # histograms are identical
identical(a1, a2[match(a1, a2)]) # vectors are identical after matching
#> [1] TRUE
# compute accumulation rates for an age where there are NAs in MCMC samples
a1 <- rbacon::accrate.age(4500)
a2 <- accrate.age_new(4500)
plot(a1 ~ a2[a2 %in% a2[match(a1, a2)]]) # would expect all points on the identity line
hist(a1) # histograms are identical
identical(a1, a2[match(a1, a2)]) # vectors are identical after matching
#> [1] TRUE
Created on 2022-02-09 by the reprex package (v2.0.1)