v1.5.0 adding jackknife method

Project.toml

@@ -1,7 +1,7 @@
 name = "JDQMCMeasurements"
 uuid = "387dbc25-93e3-4e71-ad28-7b91e77c6da4"
 authors = ["Benjamin Cohen-Stead <[email protected]>"]
-version = "1.4.0"
+version = "1.5.0"
 
 [deps]
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
@@ -10,14 +10,16 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 ShiftedArrays = "1277b4bf-5013-50f5-be3d-901d8477a67a"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 FFTW = "1.5"
 LatticeUtilities = "2.0.5"
+OffsetArrays = "1.12"
 ShiftedArrays = "2.0"
 StaticArrays = "1.5"
-OffsetArrays = "1.12"
 julia = "1.8"
+Statistics = "1.8"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

docs/src/api.md

@@ -39,11 +39,13 @@ measure_Nsqrd
 - [`fourier_transform!`](@ref)
 - [`susceptibility!`](@ref)
 - [`susceptibility`](@ref)
+- [`jackknife`](@ref)
 
 ```@docs
 fourier_transform!
 susceptibility!
 susceptibility
+jackknife
 ```
 
 ## Developer API

src/JDQMCMeasurements.jl

@@ -1,6 +1,7 @@
 module JDQMCMeasurements
 
 using LinearAlgebra
+using Statistics
 using FFTW
 using StaticArrays
 using LatticeUtilities
@@ -50,4 +51,8 @@ export bond_correlation!
 include("correlation_measurements/current_correlation.jl")
 export current_correlation!
 
+# implements the jackknife algorithm
+include("jackknife.jl")
+export jackknife
+
 end

src/jackknife.jl

@@ -0,0 +1,61 @@
+@doc raw"""
+    jackknife(
+        g::Function,
+        samples...;
+        # KEYWORD ARGUMENTS
+        bias_corrected = true,
+        jackknife_samples = similar.(samples),
+        jackknife_g = similar(samples[1])
+    )
+
+Propagate errors through the evaluation of a function `g` given the binned `samples`,
+returning both the mean and error.
+If the keyword argument `bias = true`, then the ``\mathcal{O}(1/N)`` bias is corrected.
+The keyword arguments `jackknife_samples` and `jackknife_g` can be passed to avoid
+temporary memory allocations.
+"""
+function jackknife(
+    g::Function,
+    samples...;
+    # KEYWORD ARGUMENTS
+    bias_corrected = true,
+    jackknife_samples = similar.(samples),
+    jackknife_g = similar(samples[1])
+)
+
+    # get sample size
+    N = length(jackknife_g)
+
+    # iterate over input variables
+    for i in eachindex(samples)
+
+        # calculate mean current input variable
+        x̄ = mean(samples[i])
+
+        # iterate over samples
+        for j in eachindex(samples[i])
+
+            # calculate jackknife sample by updating the mean to
+            # reflect removing the j'th sample
+            jackknife_samples[i][j] = x̄ + (x̄ - samples[i][j])/(N-1)
+        end
+    end
+
+    # evaluate the input function for each jackknife sample
+    @. jackknife_g = g(jackknife_samples...)
+
+    # calculate jackknife mean
+    ḡ = mean(jackknife_g)
+
+    # calculate jackkife error
+    Δg = sqrt( (N-1) * varm(jackknife_g, ḡ, corrected=false) )
+
+    # correct O(1/N) bias, usually doesn't matter as error scales as O(1/sqrt(N))
+    # and is typically much larger than the bias
+    if bias_corrected
+        Ḡ = g(map(mean, samples)...)
+        ḡ = N * Ḡ - (N-1) * ḡ
+    end
+
+    return ḡ, Δg
+end