Make bounding slightly more conservative given reported uncertainties

Author: brenhinkeller

src/StratMetropolis.jl

@@ -53,6 +53,8 @@
             # If bounding is requested, add extrapolated top and bottom bounds to avoid
             # issues with the stratigraphic markov chain wandering off to +/- infinity
             (youngest, oldest) = nanextrema([nanminimum(Age[Age_Sidedness .== -1]); Age[Age_Sidedness .== 0]; nanmaximum(Age[Age_Sidedness .== 1])])
+            oldest += aveuncert*1.1284/2
+            youngest -= aveuncert*1.1284/2
             dt_dH = (oldest-youngest)/(top-bottom)
             offset = round((top-bottom)*bounding/resolution)*resolution
             Age = [oldest + offset*dt_dH; Age; youngest - offset*dt_dH]
@@ -117,6 +119,8 @@
             # If bounding is requested, add extrapolated top and bottom bounds to avoid
             # issues with the stratigraphic markov chain wandering off to +/- infinity
             (youngest, oldest) = nanextrema([nanminimum(Age[Age_Sidedness .== -1]); Age[Age_Sidedness .== 0]; nanmaximum(Age[Age_Sidedness .== 1])])
+            oldest += aveuncert*1.1284/2
+            youngest -= aveuncert*1.1284/2
             dt_dH = (oldest-youngest)/(top-bottom)
             offset = round((top-bottom)*bounding/resolution)*resolution
             Age = [oldest + offset*dt_dH; Age; youngest - offset*dt_dH]
@@ -180,7 +184,9 @@
         if bounding>0
             # If bounding is requested, add extrapolated top and bottom bounds to avoid
             # issues with the stratigraphic markov chain wandering off to +/- infinity
-            (youngest, oldest) = extrema(mean.(ages))
+            (youngest, oldest) = nanextrema([nanminimum(mean.(ages[Age_Sidedness .== -1])); mean.(ages[Age_Sidedness .== 0]); nanmaximum(mean.(ages[Age_Sidedness .== 1]))])
+            oldest += aveuncert*1.1284/2
+            youngest -= aveuncert*1.1284/2
             dt_dH = (oldest-youngest)/(top-bottom)
             offset = round((top-bottom)*bounding/resolution)*resolution
             ages = unionize([Normal(oldest+offset*dt_dH, aveuncert/10); 
@@ -244,6 +250,8 @@
             # If bounding is requested, add extrapolated top and bottom bounds to avoid
             # issues with the stratigraphic markov chain wandering off to +/- infinity
             (youngest, oldest) = extrema(mean.(ages))
+            oldest += aveuncert*1.1284/2
+            youngest -= aveuncert*1.1284/2
             dt_dH = (oldest-youngest)/(top-bottom)
             offset = round((top-bottom)*bounding/resolution)*resolution
             ages = unionize([Normal(oldest+offset*dt_dH, aveuncert/10); 
@@ -333,6 +341,8 @@
             # If bounding is requested, add extrapolated top and bottom bounds to avoid
             # issues with the stratigraphic markov chain wandering off to +/- infinity
             (youngest, oldest) = nanextrema([nanminimum(Age[Age_Sidedness .== -1]); Age[Age_Sidedness .== 0]; nanmaximum(Age[Age_Sidedness .== 1])])
+            oldest += aveuncert*1.1284/2
+            youngest -= aveuncert*1.1284/2
             dt_dH = (oldest-youngest)/(top-bottom)
             offset = round((top-bottom)*bounding/resolution)*resolution
             Age = [oldest + offset*dt_dH; Age; youngest - offset*dt_dH]
@@ -401,6 +411,8 @@
             # If bounding is requested, add extrapolated top and bottom bounds to avoid
             # issues with the stratigraphic markov chain wandering off to +/- infinity
             (youngest, oldest) = nanextrema([nanminimum(Age[Age_Sidedness .== -1]); Age[Age_Sidedness .== 0]; nanmaximum(Age[Age_Sidedness .== 1])])
+            oldest += aveuncert*1.1284/2
+            youngest -= aveuncert*1.1284/2
             dt_dH = (oldest-youngest)/(top-bottom)
             offset = round((top-bottom)*bounding/resolution)*resolution
             Age = [oldest + offset*dt_dH; Age; youngest - offset*dt_dH]
@@ -492,6 +504,8 @@
             # If bounding is requested, add extrapolated top and bottom bounds to avoid
             # issues with the stratigraphic markov chain wandering off to +/- infinity
             (youngest, oldest) = nanextrema([nanminimum(Age[Age_Sidedness .== -1]); Age[Age_Sidedness .== 0]; nanmaximum(Age[Age_Sidedness .== 1])])
+            oldest += aveuncert*1.1284/2
+            youngest -= aveuncert*1.1284/2
             dt_dH = (oldest-youngest)/(top-bottom)
             offset = round((top-bottom)*bounding/resolution)*resolution
             Age = [oldest + offset*dt_dH; Age; youngest - offset*dt_dH]
@@ -560,6 +574,8 @@
             # If bounding is requested, add extrapolated top and bottom bounds to avoid
             # issues with the stratigraphic markov chain wandering off to +/- infinity
             (youngest, oldest) = nanextrema([nanminimum(Age[Age_Sidedness .== -1]); Age[Age_Sidedness .== 0]; nanmaximum(Age[Age_Sidedness .== 1])])
+            oldest += aveuncert*1.1284/2
+            youngest -= aveuncert*1.1284/2
             dt_dH = (oldest-youngest)/(top-bottom)
             offset = round((top-bottom)*bounding/resolution)*resolution
             Age = [oldest + offset*dt_dH; Age; youngest - offset*dt_dH]

test/testCoupled.jl

@@ -45,8 +45,8 @@ println("StratMetropolisDist:")
 # Test that results match expectation, within some tolerance
 @test mdl.Age isa Vector{Float64}
 @test mdl.Age ≈ [66.07, 66.06, 66.05, 66.03, 66.02, 66.01, 65.99, 65.98, 65.97, 65.96, 65.94, 65.94, 65.93, 65.93, 65.89] atol=0.1
-@test mdl.Age_025CI ≈ [66.04, 66.0, 65.98, 65.96, 65.95, 65.94, 65.94, 65.93, 65.93, 65.92, 65.92, 65.91, 65.9, 65.89, 65.82] atol=0.15
-@test mdl.Age_975CI ≈ [66.1, 66.09, 66.09, 66.08, 66.08, 66.07, 66.06, 66.05, 66.03, 66.01, 65.97, 65.97, 65.96, 65.96, 65.95] atol=0.15
+@test mdl.Age_025CI ≈ [66.04, 66.0, 65.98, 65.96, 65.95, 65.94, 65.94, 65.93, 65.93, 65.92, 65.92, 65.91, 65.9, 65.89, 65.82] atol=0.3
+@test mdl.Age_975CI ≈ [66.1, 66.09, 66.09, 66.08, 66.08, 66.07, 66.06, 66.05, 66.03, 66.01, 65.97, 65.97, 65.96, 65.96, 65.95] atol=0.3
 # Test that all age-depth models are in stratigraphic order
 @test all([issorted(x, rev=true) for x in eachcol(agedist)])
 @test all(!isnan, agedist)
@@ -66,8 +66,8 @@ println("StratMetropolisDist with systematic uncertainties:")
 # Test that results match expectation, within some tolerance
 @test mdl.Age isa Vector{Float64}
 @test mdl.Age ≈ [66.07, 66.06, 66.05, 66.03, 66.02, 66.01, 65.99, 65.98, 65.97, 65.95, 65.94, 65.94, 65.93, 65.93, 65.89] atol=0.1
-@test mdl.Age_025CI ≈ [66.04, 66.0, 65.98, 65.96, 65.95, 65.94, 65.93, 65.93, 65.92, 65.92, 65.92, 65.91, 65.9, 65.88, 65.82] atol=0.15
-@test mdl.Age_975CI ≈ [66.1, 66.09, 66.09, 66.08, 66.08, 66.07, 66.06, 66.05, 66.03, 66.01, 65.97, 65.97, 65.96, 65.96, 65.95] atol=0.15
+@test mdl.Age_025CI ≈ [66.04, 66.0, 65.98, 65.96, 65.95, 65.94, 65.93, 65.93, 65.92, 65.92, 65.92, 65.91, 65.9, 65.88, 65.82] atol=0.3
+@test mdl.Age_975CI ≈ [66.1, 66.09, 66.09, 66.08, 66.08, 66.07, 66.06, 66.05, 66.03, 66.01, 65.97, 65.97, 65.96, 65.96, 65.95] atol=0.3
 # Test that all age-depth models are in stratigraphic order
 @test all([issorted(x, rev=true) for x in eachcol(agedist)])
 @test all(!isnan, agedist)
@@ -90,8 +90,8 @@ println("StratMetropolisDist with hiata:")
 # Test that results match expectation, within some tolerance
 @test mdl.Age isa Vector{Float64}
 @test mdl.Age ≈ [66.08, 66.08, 66.07, 66.07, 66.07, 66.02, 66.02, 66.01, 66.01, 65.94, 65.94, 65.93, 65.93, 65.92, 65.89] atol=0.1
-@test mdl.Age_025CI ≈ [66.06, 66.05, 66.04, 66.03, 66.02, 65.94, 65.94, 65.94, 65.93, 65.91, 65.91, 65.9, 65.89, 65.88, 65.82] atol=0.15
-@test mdl.Age_975CI ≈ [66.11, 66.1, 66.1, 66.1, 66.1, 66.08, 66.08, 66.08, 66.08, 65.98, 65.96, 65.96, 65.96, 65.96, 65.95] atol=0.15
+@test mdl.Age_025CI ≈ [66.06, 66.05, 66.04, 66.03, 66.02, 65.94, 65.94, 65.94, 65.93, 65.91, 65.91, 65.9, 65.89, 65.88, 65.82] atol=0.3
+@test mdl.Age_975CI ≈ [66.11, 66.1, 66.1, 66.1, 66.1, 66.08, 66.08, 66.08, 66.08, 65.98, 65.96, 65.96, 65.96, 65.96, 65.95] atol=0.3
 # Test that all age-depth models are in stratigraphic order
 @test all([issorted(x, rev=true) for x in eachcol(agedist)])
 @test all(!isnan, agedist)
@@ -114,8 +114,8 @@ println("StratMetropolisDist with fitted Gaussians:")
 @time (mdl, agedist, lldist) = StratMetropolisDist(smpl, config)
 @test mdl.Age isa Vector{Float64}
 @test mdl.Age ≈ [65.99, 65.98, 65.97, 65.96, 65.96, 65.95, 65.94, 65.93, 65.92, 65.91, 65.91, 65.9, 65.88, 65.86, 65.84] atol=0.1
-@test mdl.Age_025CI ≈ [65.87, 65.86, 65.86, 65.85, 65.85, 65.84, 65.84, 65.83, 65.83, 65.83, 65.82, 65.8, 65.75, 65.73, 65.7] atol=0.15
-@test mdl.Age_975CI ≈ [66.09, 66.09, 66.08, 66.08, 66.08, 66.07, 66.06, 66.05, 66.04, 66.02, 66.0, 65.99, 65.98, 65.97, 65.96] atol=0.15
+@test mdl.Age_025CI ≈ [65.87, 65.86, 65.86, 65.85, 65.85, 65.84, 65.84, 65.83, 65.83, 65.83, 65.82, 65.8, 65.75, 65.73, 65.7] atol=0.3
+@test mdl.Age_975CI ≈ [66.09, 66.09, 66.08, 66.08, 66.08, 66.07, 66.06, 66.05, 66.04, 66.02, 66.0, 65.99, 65.98, 65.97, 65.96] atol=0.3
 # Test that all age-depth models are in stratigraphic order
 @test all([issorted(x, rev=true) for x in eachcol(agedist)])
 @test all(!isnan, agedist)

test/testStratOnly.jl

@@ -19,7 +19,7 @@ config.resolution = 5 # Same units as sample height. Smaller is slower!
 config.bounding = 0.5 # how far away do we place runaway bounds, as a fraction of total section height. Larger is slower.
 (bottom, top) = extrema(smpl.Height)
 npoints_approx = round(Int,length(bottom:config.resolution:top) * (1 + 2*config.bounding))
-config.nsteps = 100000 # Number of steps to run in distribution MCMC
+config.nsteps = 1000000 # Number of steps to run in distribution MCMC
 config.burnin = 100000*npoints_approx # Number to discard
 config.sieve = round(Int,npoints_approx) # Record one out of every nsieve steps
 
@@ -28,9 +28,9 @@ config.sieve = round(Int,npoints_approx) # Record one out of every nsieve steps
 
 # Test that results match expectation, within some tolerance
 @test mdl.Age isa Vector{Float64}
-@test mdl.Age ≈  [774.446455402488, 764.138871634567, 753.6234690338521, 743.7736679540459, 734.1288211670727, 724.4222568319159, 719.9641790169545, 715.7853469234702, 711.659901307201, 709.0927902799843, 706.6073112860937, 704.1349331543528, 701.576806557881, 698.7594637725156, 695.9077192905107] atol=1
-@test mdl.Age_025CI ≈ [749.7477463948217, 746.630273097855, 743.9690750302638, 723.4084389968576, 718.3588586689344, 715.5410003896017, 706.7023568283762, 703.1283100699234, 700.9213523438647, 697.904849247825, 696.2105445428034, 695.0433927602111, 694.1016935971145, 693.2625238063798, 684.8895558316834] atol=3
-@test mdl.Age_975CI ≈ [806.7021880750214, 794.9197023904139, 763.2743319392429, 760.0094083209071, 754.7169964175238, 733.5482560401001, 731.2806438300145, 728.3050781269916, 722.3491038457139, 720.830748668773, 718.9846443130959, 716.6641608791327, 713.2397206761551, 704.3594827974073, 703.4763901911105] atol=3
+@test mdl.Age ≈  [775.2009717531292, 764.5613873063818, 753.7193603935707, 743.859942276044, 734.1903641532622, 724.4313056816793, 719.9612277736258, 715.7765014002005, 711.6562087664921, 709.0844186854755, 706.5904994086322, 704.1094515386092, 701.5311095414545, 698.7091903065589, 695.7689678220307]  atol=1
+@test mdl.Age_025CI ≈ [750.0939493479998, 746.7621661695302, 744.0833309387968, 723.3669863534756, 718.3636351303046, 715.540624074875, 706.7137614573669, 703.1038143674268, 700.844732326315, 697.8358269342763, 696.1454906755058, 694.9648114134607, 694.0249776716082, 693.1542101854486, 684.3681541813636] atol=3
+@test mdl.Age_975CI ≈ [807.2510236903488, 795.3902664676122, 763.3234305018414, 760.0376619032745, 754.7160594008493, 733.5082697327649, 731.2581768641392, 728.2692821103202, 722.3718766353179, 720.8356091969036, 718.9801564927342, 716.6713426324544, 713.1583260608256, 704.316125432296, 703.4023532440506] atol=3
 # Test that all age-depth models are in stratigraphic order
 @test all([issorted(x, rev=true) for x in eachcol(agedist)])
 @test all(!isnan, agedist)
@@ -51,12 +51,12 @@ hiatus.Duration_sigma = [   3.1,    2.0 ]
 
 # Test that results match expectation, within some tolerance
 @test mdl.Age isa Vector{Float64}
-@test mdl.Age ≈ [774.9621590393864, 764.7483478838244, 754.3928287806855, 749.1279876514806, 729.1819943547974, 724.1793043771663, 720.8913053718023, 717.7627759721249, 714.5509500060313, 713.2002472528285, 702.1296086101366, 700.7704707147774, 699.4643412956174, 698.0733121394825, 695.3696935693185] atol=1
-@test mdl.Age_025CI ≈ [750.6328462444853, 747.6305055671014, 745.2219453018333, 735.8795895625685, 717.8606462193696, 715.6849749469166, 710.1425439567424, 707.5591665793866, 705.808337538727, 704.4299789235162, 694.3410849948168, 693.699614909472, 693.1108571564564, 692.5485716537205, 684.6340487749085] atol=3
-@test mdl.Age_975CI ≈ [806.6715104282963, 794.8753648212958, 763.7844862197752, 761.1148678026058, 742.5015528417246, 732.771524797932, 730.9079222710706, 728.4137307378472, 723.6929471023955, 722.6480626953821, 711.3953411778708, 709.6936433474348, 707.3968615013218, 703.6205967015544, 702.7858695444451] atol=3
+@test mdl.Age ≈   [775.5494538893379, 765.0398661222231, 754.4553671556596, 749.1165989011665, 729.1802542341943, 724.1690788281638, 720.8671584198628, 717.7310672911206, 714.5176266505601, 713.1589898075201, 702.0799302600404, 700.7161043228054, 699.4056554863363, 698.0132290102006, 695.2079260036835] atol=1
+@test mdl.Age_025CI ≈  [750.8736692043892, 747.7489306192067, 745.269061071627, 735.8291026862441, 717.8684206092925, 715.6903606224124, 710.0972278697697, 707.5386715468874, 705.8103447794734, 704.3928559634423, 694.2857415733332, 693.6245799703629, 693.0486461983832, 692.4919675548376, 684.2374067953986] atol=3
+@test mdl.Age_975CI ≈  [807.3229235373376, 795.5005210793569, 763.8089856382381, 761.1489325004852, 742.73112549815, 732.7647859218115, 730.8975324858483, 728.4214438533553, 723.6903413737293, 722.624657136144, 711.3505494727046, 709.6442679173194, 707.3487791460136, 703.5758337472577, 702.7247254673973] atol=3
 # Test that all age-depth models are in stratigraphic order
 @test all([issorted(x, rev=true) for x in eachcol(agedist)])
 @test all(!isnan, agedist)
 @test size(hiatusdist) == (nhiatuses, config.nsteps)
-@test mean(hiatusdist, dims=2) ≈ [10.580012942504894; 18.96167245288326;;] atol=2
+@test mean(hiatusdist, dims=2) ≈ [11.083103050893387; 19.955414496172235;;] atol=2
 @test -Inf < mean(lldist) < 0