Increase default cutoff in distribution bootstrapping given better empirical performance

Author: brenhinkeller

Project.toml

@@ -1,7 +1,7 @@
 name = "Chron"
 uuid = "68885b1f-77b5-52a7-b2e7-6a8014c56b98"
 authors = ["C. Brenhin Keller <[email protected]>"]
-version = "0.6.5"
+version = "0.7.0"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

src/DistMetropolis.jl

@@ -2,7 +2,7 @@
 
     """
     ```julia
-    BootstrapCrystDistributionKDE(smpl::ChronAgeData; cutoff=-0.05, [tpbloss=0])
+    BootstrapCrystDistributionKDE(smpl::ChronAgeData; cutoff=0.02, [tpbloss=0])
     ```
     Bootstrap an estimate of the pre-eruptive (or pre-depositional) mineral
     crystallization distribution shape from a Chron.ChronAgeData object containing
@@ -20,7 +20,7 @@
     BootstrappedDistribution = BootstrapCrystDistributionKDE(smpl)
     ```
     """
-    function BootstrapCrystDistributionKDE(smpl::ChronAgeData; cutoff::Number=-0.05, tpbloss::Number=0)
+    function BootstrapCrystDistributionKDE(smpl::ChronAgeData; cutoff::Number=0.02, tpbloss::Number=0)
         # Extact variables froms struct
         Name = collect(smpl.Name)::Vector{String}
         Path = smpl.Path::String
@@ -66,7 +66,7 @@
             end
         end
 
-        # Calculate kernel density estimate, truncated at 0
+        # Calculate kernel density estimate, truncated at `cutoff`
         kd = kde(allscaled,npoints=2^7)
         t = kd.x .> cutoff # Ensure sharp cutoff at eruption / deposition
         return kd.density[t]
@@ -75,7 +75,7 @@
 
     """
     ```julia
-    BootstrapCrystDistributionKDE(data::AbstractArray, [sigma::AbstractArray]; cutoff=-0.05)
+    BootstrapCrystDistributionKDE(data::AbstractArray, [sigma::AbstractArray]; cutoff=0.02)
     ```
     Bootstrap an estimate of the pre-eruptive (or pre-depositional) mineral
     crystallization distribution shape from a 1- or 2-d array of sample ages
@@ -90,7 +90,7 @@
     BootstrappedDistribution = BootstrapCrystDistributionKDE(1:10, ones(10))
     ```
     """
-    function BootstrapCrystDistributionKDE(data::AbstractArray{T}; cutoff::Number=-0.05) where {T<:Number}
+    function BootstrapCrystDistributionKDE(data::AbstractArray{T}; cutoff::Number=0.02) where {T<:Number}
         # Load all data points and scale from 0 to 1
         allscaled = Vector{float(T)}()
         for i=1:size(data,2)
@@ -101,13 +101,13 @@
             append!(allscaled, scaled)
         end
 
-        # Calculate kernel density estimate, truncated at 0
+        # Calculate kernel density estimate, truncated at `cutoff`
         kd = kde(allscaled,npoints=2^7)
         t = kd.x .> cutoff
         return kd.density[t]
     end
 
-    function BootstrapCrystDistributionKDE(data::AbstractArray{T}, sigma::AbstractArray{<:Number}; cutoff::Number=-0.05) where {T<:Number}
+    function BootstrapCrystDistributionKDE(data::AbstractArray{T}, sigma::AbstractArray{<:Number}; cutoff::Number=0.02) where {T<:Number}
         # Array to hold stacked, scaled data
         allscaled = Vector{float(T)}()
 
@@ -131,7 +131,7 @@
             end
         end
 
-        # Calculate kernel density estimate, truncated at 0
+        # Calculate kernel density estimate, truncated at `cutoff`
         kd = kde(allscaled,npoints=2^7)
         t = kd.x .> cutoff
         return kd.density[t]

test/testCoupled.jl

@@ -17,7 +17,7 @@ smpl = screen_outliers(smpl, maxgap=50; make_plots)
 # Distribution boostrapping from chron strat object
 BootstrappedDistribution = BootstrapCrystDistributionKDE(smpl)
 @test BootstrappedDistribution isa Vector{Float64}
-@test BootstrappedDistribution[1] ≈ 1.0511396290122654
+@test BootstrappedDistribution[1] ≈ 1.2986370127270708
 
 # Estimate the eruption age distributions for each sample  - - - - - - - -
 
@@ -138,7 +138,7 @@ smpl = screen_outliers(smpl, maxgap=100; make_plots, discordancemin=-5)
 
 BootstrappedDistribution = BootstrapCrystDistributionKDE(smpl, tpbloss=0)
 @test BootstrappedDistribution isa Vector{Float64}
-@test BootstrappedDistribution[1] ≈ 0.3697041339884247
+@test BootstrappedDistribution[1] ≈ 0.45349905011835107
 
 # Configure distribution model here
 distSteps = 2*10^5 # Number of steps to run in distribution MCMC

test/testDist.jl

@@ -4,22 +4,30 @@
     dist = MeltsVolcanicZirconDistribution ./ nanmean(MeltsVolcanicZirconDistribution)
     @test Chron.Isoplot.dist_ll(dist, mu, sigma, 100,101) ≈ -3.6933357793356576
 
-    tmindist = metropolis_min(2*10^5, dist, mu, sigma, burnin=10^5)
+    tmindist = metropolis_min(2*10^5, MeltsVolcanicZirconDistribution, mu, sigma, burnin=10^5)
     @test isapprox(nanmean(tmindist), 99.925, atol=0.015)
-
+    
     tmindist, tmaxdist, lldist, acceptancedist = metropolis_minmax(2*10^5, MeltsVolcanicZirconDistribution, mu, sigma, burnin=10^5)
     @test isapprox(nanmean(tmindist), 99.925, atol=0.015)
     @test isapprox(nanmean(tmaxdist), 101.08, atol=0.015)
     @test -Inf < mean(lldist) < 0
     @test 0 < mean(acceptancedist) < 1
 
+    tmindist = metropolis_min(2*10^5, UniformDistribution, mu, sigma, burnin=10^5)
+    @test isapprox(nanmean(tmindist), 99.961, atol=0.015)
+
 ## -- Test distribution bootstrapping
 
     dist = BootstrapCrystDistributionKDE(1:10)
-    @test all(isapprox.(dist, [0.465404, 0.503149, 0.540308, 0.576485, 0.611314, 0.644472, 0.675691, 0.704756, 0.731516, 0.75588, 0.777815, 0.797345, 0.81454, 0.82951, 0.842398, 0.85337, 0.862605, 0.87029, 0.876611, 0.881749, 0.885873, 0.889136, 0.891676, 0.893612, 0.895043, 0.896047, 0.896684, 0.896993, 0.896993, 0.896684, 0.896047, 0.895043, 0.893612, 0.891676, 0.889136, 0.885873, 0.881749, 0.876611, 0.87029, 0.862605, 0.85337, 0.842398, 0.82951, 0.81454, 0.797345, 0.777815, 0.75588, 0.731516, 0.704756, 0.675691, 0.644472, 0.611314, 0.576485, 0.540308, 0.503149, 0.465404, 0.427489, 0.389826, 0.352829, 0.316888, 0.282361, 0.249559, 0.218742, 0.19011, 0.163804, 0.139901, 0.118424, 0.099338, 0.082567, 0.067992, 0.055465, 0.044819, 0.03587, 0.028432, 0.022318, 0.017347, 0.01335, 0.010173, 0.007675, 0.005732, 0.004238, 0.003102, 0.002247, 0.001612, 0.001145, 0.000806, 0.000563, 0.000392, 0.000275, 0.000197, 0.00015, 0.000128], atol=0.00001))
+    @test all(isapprox.(dist, [0.611314, 0.644472, 0.675691, 0.704756, 0.731516, 0.75588, 0.777815, 0.797345, 0.81454, 0.82951, 0.842398, 0.85337, 0.862605, 0.87029, 0.876611, 0.881749, 0.885873, 0.889136, 0.891676, 0.893612, 0.895043, 0.896047, 0.896684, 0.896993, 0.896993, 0.896684, 0.896047, 0.895043, 0.893612, 0.891676, 0.889136, 0.885873, 0.881749, 0.876611, 0.87029, 0.862605, 0.85337, 0.842398, 0.82951, 0.81454, 0.797345, 0.777815, 0.75588, 0.731516, 0.704756, 0.675691, 0.644472, 0.611314, 0.576485, 0.540308, 0.503149, 0.465404, 0.427489, 0.389826, 0.352829, 0.316888, 0.282361, 0.249559, 0.218742, 0.19011, 0.163804, 0.139901, 0.118424, 0.099338, 0.082567, 0.067992, 0.055465, 0.044819, 0.03587, 0.028432, 0.022318, 0.017347, 0.01335, 0.010173, 0.007675, 0.005732, 0.004238, 0.003102, 0.002247, 0.001612, 0.001145, 0.000806, 0.000563, 0.000392, 0.000275, 0.000197, 0.00015, 0.000128], atol=0.00001))
 
     dist = BootstrapCrystDistributionKDE(1:10, ones(10))
-    @test all(isapprox.(dist, [0.500228, 0.535143, 0.569056, 0.601647, 0.632641, 0.661805, 0.688961, 0.713982, 0.736794, 0.757375, 0.775748, 0.79198, 0.806169, 0.818443, 0.828947, 0.837842, 0.845291, 0.851459, 0.856502, 0.86057, 0.863796, 0.866299, 0.868179, 0.869517, 0.870375, 0.870794, 0.870794, 0.870375, 0.869517, 0.868179, 0.866299, 0.863796, 0.86057, 0.856502, 0.851459, 0.845291, 0.837842, 0.828947, 0.818443, 0.806169, 0.79198, 0.775748, 0.757375, 0.736794, 0.713982, 0.688961, 0.661805, 0.632641, 0.601647, 0.569056, 0.535143, 0.500228, 0.464656, 0.428798, 0.39303, 0.357728, 0.323251, 0.289935, 0.258079, 0.227938, 0.199721, 0.173583, 0.149624, 0.127894, 0.108393, 0.091075, 0.075858, 0.062627, 0.051244, 0.041553, 0.03339, 0.026585, 0.020972, 0.016391, 0.012691, 0.009734, 0.007396, 0.005566, 0.004149, 0.003064, 0.002241, 0.001623, 0.001166, 0.00083, 0.000587, 0.000415, 0.000295, 0.000215, 0.000166, 0.000143], atol=0.00001))
+    @test all(isapprox.(dist, [0.601647, 0.632641, 0.661805, 0.688961, 0.713982, 0.736794, 0.757375, 0.775748, 0.79198, 0.806169, 0.818443, 0.828947, 0.837842, 0.845291, 0.851459, 0.856502, 0.86057, 0.863796, 0.866299, 0.868179, 0.869517, 0.870375, 0.870794, 0.870794, 0.870375, 0.869517, 0.868179, 0.866299, 0.863796, 0.86057, 0.856502, 0.851459, 0.845291, 0.837842, 0.828947, 0.818443, 0.806169, 0.79198, 0.775748, 0.757375, 0.736794, 0.713982, 0.688961, 0.661805, 0.632641, 0.601647, 0.569056, 0.535143, 0.500228, 0.464656, 0.428798, 0.39303, 0.357728, 0.323251, 0.289935, 0.258079, 0.227938, 0.199721, 0.173583, 0.149624, 0.127894, 0.108393, 0.091075, 0.075858, 0.062627, 0.051244, 0.041553, 0.03339, 0.026585, 0.020972, 0.016391, 0.012691, 0.009734, 0.007396, 0.005566, 0.004149, 0.003064, 0.002241, 0.001623, 0.001166, 0.00083, 0.000587, 0.000415, 0.000295, 0.000215, 0.000166, 0.000143], atol=0.00001))
+
+    mu, sigma = collect(100:0.1:101), 0.01*ones(11);
+    tmindist = metropolis_min(2*10^5, dist, mu, sigma, burnin=10^5)
+    @test isapprox(nanmean(tmindist), 99.961, atol=0.015)
+
 
 ## -- Test systematic uncertainty propagation