Use histmean and histstd to determine parameters of Radiocarbon…

… distributions

Project.toml

@@ -22,7 +22,7 @@ Distributions = "0.15 - 0.25"
 Isoplot = "0.3"
 KernelDensity = "0.4 - 0.6"
 LsqFit = "0.7 - 0.13"
-NaNStatistics = "0.4, 0.5, 0.6"
+NaNStatistics = "0.6.38"
 Plots = "1"
 ProgressMeter = "1"
 Reexport = "0.2, 1.0"

src/Utilities.jl

@@ -107,10 +107,10 @@
 
         dist = exp.(ldist)
         dist ./= sum(dist) * calibration.dt # Normalize
-        samples = draw_from_distribution(dist, 10^6) .* maximum(calibration.Age_Calendar)
-        μ, σ = nanmean(samples), nanstd(samples)
+        μ = histmean(dist, calibration.Age_Calendar)
+        σ = histstd(dist, calibration.Age_Calendar, corrected=false)
 
-        return Radiocarbon(μ, σ, ldist)
+        return Radiocarbon(μ, σ, dist, ldist)
     end
 
     ## Conversions
@@ -135,6 +135,10 @@
         return linterp_at_index(d.ldist, x, -maxintfloat(T))
     end
 
+    ## Statistics
+    Distributions.mean(d::Radiocarbon) = d.μ
+    Distributions.std(d::Radiocarbon) = d.σ
+
 ## --- 
 
     # Interpolate log likelihood from an array

test/testUtilities.jl

@@ -24,8 +24,11 @@
     @test pdf.(d, [900, 950, 1000]) ≈ [3.372067127114498e-7, 0.0036210273644940918, 1.5582187464154278e-12]
 
     @test eltype(d) === partype(d) === Float64
-    @test location(d) ≈ 926.3094740428785 atol=0.1
-    @test scale(d) ≈ 7.411299532288234 atol=0.1
+    @test location(d) ≈ 927.2875322569857
+    @test scale(d) ≈ 7.412902965559571
+
+    @test mean(d) ≈ 927.2875322569857
+    @test std(d) ≈ 7.412902965559571
 
 ## --- Other utility functions for log likelihoods