diff --git a/EC.jl b/EC.jl
index 61cadb6..adfd866 100644
--- a/EC.jl
+++ b/EC.jl
@@ -1,4 +1,4 @@
-using SparseArrays, LinearAlgebra, Arpack, Plots
+using Statistics, SparseArrays, LinearAlgebra, Arpack, Plots
 include("common.jl")
 
 "EC model for a Hamiltonian family H(c) = H0 + c * H1"
@@ -12,11 +12,18 @@ mutable struct affine_EC
     H1_EC
     N_EC
 
+    ensemble_size::Int
+    H0_EC_ensemble
+    H1_EC_ensemble
+    N_EC_ensemble
+
     training_E::Vector{ComplexF64}
     exact_E::Vector{ComplexF64}
     extrapolated_E::Vector{ComplexF64}
+    extrapolated_CI::Vector{ComplexF64}
 
-    affine_EC(H0::AbstractMatrix{ComplexF64}, H1::AbstractMatrix{ComplexF64}, weights::Vector{ComplexF64}=ones(ComplexF64, size(H0, 1))) = new(H0, H1, weights, false, nothing, nothing, nothing, ComplexF64[], ComplexF64[], ComplexF64[])
+    affine_EC(H0::AbstractMatrix{ComplexF64}, H1::AbstractMatrix{ComplexF64}, weights::Vector{ComplexF64}=ones(ComplexF64, size(H0, 1)); ensemble_size=0) =
+        new(H0, H1, weights, false, nothing, nothing, nothing, ensemble_size, Matrix[], Matrix[], Matrix[], ComplexF64[], ComplexF64[], ComplexF64[], ComplexF64[])
 end
 
 "Train an EC model for a given range of c values.
@@ -51,18 +58,20 @@ function train!(EC::affine_EC, c_vals; ref_eval=-10.0, pseudo_inv_rtol=1e-6, gra
     end
 
     CAEC && append!(training_vecs, conj.(training_vecs))
+    weights_mat = spdiagm(EC.weights)
 
-    if gram_schmidt_threshold ≠ 0
+    if gram_schmidt_threshold > 0
+        EC.ensemble_size > 0 && generate_resampled_EC_matrices(EC, training_vecs, weights_mat, gram_schmidt_threshold)
         training_vecs = gram_schmidt!(training_vecs, EC.weights, gram_schmidt_threshold; verbose=verbose)
     end
 
     EC_basis = hcat(training_vecs...)
-    weights_mat = spdiagm(EC.weights)
     EC.H0_EC = transpose(EC_basis) * weights_mat * EC.H0 * EC_basis
     EC.H1_EC = transpose(EC_basis) * weights_mat * EC.H1 * EC_basis
     EC.N_EC = transpose(EC_basis) * weights_mat * EC_basis
 
     if pseudo_inv_rtol > 0
+        EC.ensemble_size > 0 && generate_resampled_EC_matrices(EC, pseudo_inv_rtol)
         inv_N_EC = pinv(EC.N_EC; rtol=pseudo_inv_rtol)
         EC.H0_EC = inv_N_EC * EC.H0_EC
         EC.H1_EC = inv_N_EC * EC.H1_EC
@@ -72,6 +81,38 @@ function train!(EC::affine_EC, c_vals; ref_eval=-10.0, pseudo_inv_rtol=1e-6, gra
     EC.trained = true
 end
 
+resample(n::Int) = rand(1:n, n) |> unique |> sort
+
+"Generate an resampled ensemble of reduced matrices performing Gram-Schmidt orthonormalization for each."
+function generate_resampled_EC_matrices(EC::affine_EC, training_vecs, weights_mat, gram_schmidt_threshold)
+    for _ in 1:EC.ensemble_size
+        sampled_vecs = deepcopy(training_vecs[resample(length(training_vecs))])
+        orth_vecs = gram_schmidt!(sampled_vecs, EC.weights, gram_schmidt_threshold; verbose=false)
+        EC_basis = hcat(orth_vecs...)
+
+        H0_EC = transpose(EC_basis) * weights_mat * EC.H0 * EC_basis
+        H1_EC = transpose(EC_basis) * weights_mat * EC.H1 * EC_basis
+        N_EC = transpose(EC_basis) * weights_mat * EC_basis
+        push!(EC.H0_EC_ensemble, H0_EC)
+        push!(EC.H1_EC_ensemble, H1_EC)
+        push!(EC.N_EC_ensemble, N_EC)
+    end
+end
+
+"Generate an resampled ensemble of reduced matrices performing Moore-Penrose psuedoinverse for each."
+function generate_resampled_EC_matrices(EC::affine_EC, pseudo_inv_rtol)
+    for _ in 1:EC.ensemble_size
+        sample = resample(size(EC.N_EC, 1))
+        new_N_EC = EC.N_EC[sample, sample]
+        new_H0_EC = EC.H0_EC[sample, sample]
+        new_H1_EC = EC.H1_EC[sample, sample]
+
+        inv_N_EC = pinv(new_N_EC; rtol=pseudo_inv_rtol)
+        push!(EC.H0_EC_ensemble, inv_N_EC * new_H0_EC)
+        push!(EC.H1_EC_ensemble, inv_N_EC * new_H1_EC)
+    end
+end
+
 "Extrapolate using a trained EC model for a given range of c values
 If a list is provided for ref_eval, they are used as reference values for picking the closest eigenvalues at each point.
 If a single number is provided for ref_eval, it is used as a reference for the first point, and the previous eigenvalue is used as the reference for each successive point.
@@ -106,6 +147,18 @@ function extrapolate!(EC::affine_EC, c_vals; ref_eval=EC.training_E[end], verbos
         H_EC = EC.H0_EC + c .* EC.H1_EC
         evals = isnothing(EC.N_EC) ? eigvals(H_EC) : eigvals(H_EC, EC.N_EC)
         push!(EC.extrapolated_E, nearest(evals, current_E))
+
+        if EC.ensemble_size > 0
+            E_ensemble = ComplexF64[]
+            for i in 1:EC.ensemble_size
+                H_EC = EC.H0_EC_ensemble[i] + c .* EC.H1_EC_ensemble[i]
+                evals = isempty(EC.N_EC_ensemble) ? eigvals(H_EC) : eigvals(H_EC, EC.N_EC_ensemble[i])
+                push!(E_ensemble, nearest(evals, current_E))
+            end
+            re_CI = std(real.(E_ensemble))
+            im_CI = std(imag.(E_ensemble))
+            push!(EC.extrapolated_CI, complex(re_CI, im_CI))
+        end
     end
 end
 
@@ -116,7 +169,11 @@ exportCSV(EC::affine_EC, filename) = exportCSV(filename, (EC.training_E, EC.exac
 function plot(EC::affine_EC, save_fig_filename=nothing; basis_points=nothing, basis_contour=nothing, xlims=nothing, ylims=nothing)
     scatter(real.(EC.training_E), imag.(EC.training_E), label="training")
     scatter!(real.(EC.exact_E), imag.(EC.exact_E), label="exact")
-    scatter!(real.(EC.extrapolated_E), imag.(EC.extrapolated_E), label="extrapolated")
+    if EC.ensemble_size > 0
+        scatter!(real.(EC.extrapolated_E), imag.(EC.extrapolated_E), xerror=real.(EC.extrapolated_CI), yerror=imag.(EC.extrapolated_CI), label="extrapolated")
+    else
+        scatter!(real.(EC.extrapolated_E), imag.(EC.extrapolated_E), label="extrapolated")
+    end
 
     isnothing(basis_points) || scatter!(real.(basis_points), imag.(basis_points), m=:x, label="basis")
     isnothing(basis_contour) || plot!(real.(basis_contour), imag.(basis_contour), label="contour")
diff --git a/calculations/3body_Berggren_B2R_EC.jl b/calculations/3body_Berggren_B2R_EC.jl
index 7a5f867..fbcca05 100644
--- a/calculations/3body_Berggren_B2R_EC.jl
+++ b/calculations/3body_Berggren_B2R_EC.jl
@@ -32,7 +32,7 @@ extrapolating_ref = [4.076662025307587-0.012709842443350328im,
     1.7164583929199813-0.0005455212208182736im,
     1.233088227541505-0.0003070320106485624im]
 
-EC = affine_EC(H0, Vp, weights)
+EC = affine_EC(H0, Vp, weights; ensemble_size=32)
 train!(EC, training_c; ref_eval=training_ref, CAEC=true)
 extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref)
 
diff --git a/calculations/3body_CSM_B2R_EC.jl b/calculations/3body_CSM_B2R_EC.jl
index 5de6106..51e28ee 100644
--- a/calculations/3body_CSM_B2R_EC.jl
+++ b/calculations/3body_CSM_B2R_EC.jl
@@ -32,7 +32,7 @@ exact_E = [4.076662025307587-0.012709842443350328im,
     1.7164583929199813-0.0005455212208182736im,
     1.233088227541505-0.0003070320106485624im]
 
-EC = affine_EC(H0, Vp, weights)
+EC = affine_EC(H0, Vp, weights; ensemble_size=32)
 train!(EC, training_c; ref_eval=training_ref, CAEC=true)
 extrapolate!(EC, extrapolating_c; precalculated_exact_E=exact_E)
 
diff --git a/calculations/3body_HO_B2R_EC.jl b/calculations/3body_HO_B2R_EC.jl
index 7d37c2d..7376df4 100644
--- a/calculations/3body_HO_B2R_EC.jl
+++ b/calculations/3body_HO_B2R_EC.jl
@@ -28,7 +28,7 @@ extrapolating_ref = [4.076662025307587-0.012709842443350328im,
 training_c = [2.6, 2.4, 2.2, 2.0, 1.8]
 extrapolating_c = 0.0 : 0.2 : 1.2
 
-EC = affine_EC(H0, Vp)
+EC = affine_EC(H0, Vp; ensemble_size=32)
 train!(EC, training_c; ref_eval=training_ref, CAEC=true)
 extrapolate!(EC, extrapolating_c; ref_eval=extrapolating_ref)