From bbf4334930fba36519bd748f52ef7d4852ebbca3 Mon Sep 17 00:00:00 2001
From: "C. Brenhin Keller" <cbkeller@dartmouth.edu>
Date: Thu, 27 Jun 2024 23:52:15 -0400
Subject: [PATCH] Use `histmean` and `histstd` to determine parameters of
 `Radiocarbon` distributions

---
 Project.toml          |  2 +-
 src/Utilities.jl      | 10 +++++++---
 test/testUtilities.jl |  7 +++++--
 3 files changed, 13 insertions(+), 6 deletions(-)

diff --git a/Project.toml b/Project.toml
index ad0cc44..79b1bfb 100644
--- a/Project.toml
+++ b/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"
diff --git a/src/Utilities.jl b/src/Utilities.jl
index 2825421..0852d39 100644
--- a/src/Utilities.jl
+++ b/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
diff --git a/test/testUtilities.jl b/test/testUtilities.jl
index 87e0987..f835a56 100644
--- a/test/testUtilities.jl
+++ b/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