How good are the fit uncertainty estimates?

This document looks at two different sets of measured spectra to evaluate how well, the estimated fit uncertainties match the observed uncertainties. We are using two materials known to be highly homogeneous - K412 and ADM-6006a glasses. To produce a set of spectra that differ only in count statistics, we will subdivide a spectrum into 100 spectra with an effective live-time of 0.01 of the original spectrum livetime. We expect that we should be able to compare the "fit-predicted" uncertainties with the "observed distribution" of measured values.

To be clear:

• "fit-predicted" - Comes from the covariance matrix output from the linear least squares fit
• "observed distribution" - Calculated as the standard-deviation of the 100 fit values.

We want the ratio of the (observed distribution) / (fit predicted) to be unity or close. We will call this ratio the "heterogeneity" (or "hetero" in the fourth column of the describe(....) table.)

We then repeat the process on the multiple measured unknown spectra. (4 for K412 and 15 for ADM-6005a) We expect the heterogeneity to be one or larger - larger than unity when the sample is not perfectly homogeneous.

using NeXLSpectrum              # Provides spectrum reading and fitting tools
using NeXLMatrixCorrection      # Provides `quant` to convert k-ratios to mass fraction.
using DataFrames                # Tables
using Latexify
using BenchmarkTools

K412

Load the spectra, define the fit model and apply it.

path = joinpath(@__DIR__,"K412 spectra")
unks = map(0:4) do i 
  loadspectrum(joinpath(path, "III-E K412[$i][4].msa"))
end
det = matching(unks[1], 132.0, 10)
frs = references( [
  reference(n"Al", joinpath(path, "Al2O3 std.msa"), mat"Al2O3"),
  reference(n"Mg", joinpath(path, "MgO std.msa"), mat"MgO"),
  reference(n"Fe", joinpath(path, "Fe std.msa"), mat"Fe"),
  reference(n"Si", joinpath(path, "SiO2 std.msa"), mat"SiO2"),
  reference(n"O", joinpath(path, "SiO2 std.msa"), mat"SiO2"),
  reference(n"Ca", joinpath(path, "CaF2 std.msa") ,mat"CaF2")
], det)
# frs is now a FilteredReference[] used to fit the unknowns.

# Split the counts in unks[1] into 100 randomized spectra which will sum to unks[1] then fit them
res = map(subdivide(unks[1], 100)) do s
  fit_spectrum(s, frs)
end

100-element Vector{FilterFitResult{Float64}}:
 FitResult(Sub[III-E K412[0][all],1 of 100])
 FitResult(Sub[III-E K412[0][all],2 of 100])
 FitResult(Sub[III-E K412[0][all],3 of 100])
 FitResult(Sub[III-E K412[0][all],4 of 100])
 FitResult(Sub[III-E K412[0][all],5 of 100])
 FitResult(Sub[III-E K412[0][all],6 of 100])
 FitResult(Sub[III-E K412[0][all],7 of 100])
 FitResult(Sub[III-E K412[0][all],8 of 100])
 FitResult(Sub[III-E K412[0][all],9 of 100])
 FitResult(Sub[III-E K412[0][all],10 of 100])
 ⋮
 FitResult(Sub[III-E K412[0][all],92 of 100])
 FitResult(Sub[III-E K412[0][all],93 of 100])
 FitResult(Sub[III-E K412[0][all],94 of 100])
 FitResult(Sub[III-E K412[0][all],95 of 100])
 FitResult(Sub[III-E K412[0][all],96 of 100])
 FitResult(Sub[III-E K412[0][all],97 of 100])
 FitResult(Sub[III-E K412[0][all],98 of 100])
 FitResult(Sub[III-E K412[0][all],99 of 100])
 FitResult(Sub[III-E K412[0][all],100 of 100])

Repeat the fit for the 4 measured unknowns.

res= map(unks) do s
  fit_spectrum(s, frs)
end

5-element Vector{FilterFitResult{Float64}}:
 FitResult(III-E K412[0][all])
 FitResult(III-E K412[1][all])
 FitResult(III-E K412[2][all])
 FitResult(III-E K412[3][all])
 FitResult(III-E K412[4][all])

Summary statistics.

@btime fit_spectrum(unks[1], frs)

291.700 μs (566 allocations: 261.02 KiB)
FitResult(III-E K412[0][all])

AMM-6005a

Repeat using the ADM glass. Fe is not present in ADM-6005a but we fit it to see what a null result looks like.

path = normpath(joinpath(@__DIR__, "..","test","ADM6005a spectra"))
unks = map(i->loadspectrum(joinpath(path,"ADM-6005a_$i.msa")),1:15)
al, caf2, fe, ge, si, sio2, ti, zn = map(f->loadspectrum(joinpath(path,"$f.msa")), ("Al std", "CaF2 std", "Fe std", "Ge std", "Si std", "SiO2 std", "Ti trimmed","Zn std"))

det = matching(unks[1], 132.0, 10)

frs = references( [
  reference(n"Al", al, mat"Al" ), #
  reference(n"Ca", caf2, mat"CaF2" ),   #
  reference(n"Fe", fe, mat"Fe" ),    #
  reference(n"Ge", ge, mat"Ge" ),    #
  reference(n"Si", si, mat"Si" ),  #
  reference(n"O", sio2, mat"SiO2" ),  #
  reference(n"Ti", ti, mat"Ti" ),
  reference(n"Zn", zn, mat"Zn" ) 
], det)

ss = 
res= map(subdivide(unks[1], 100)) do s
  fit_spectrum(s, frs)
end

100-element Vector{FilterFitResult{Float64}}:
 FitResult(Sub[ADM-6005a_1,1 of 100])
 FitResult(Sub[ADM-6005a_1,2 of 100])
 FitResult(Sub[ADM-6005a_1,3 of 100])
 FitResult(Sub[ADM-6005a_1,4 of 100])
 FitResult(Sub[ADM-6005a_1,5 of 100])
 FitResult(Sub[ADM-6005a_1,6 of 100])
 FitResult(Sub[ADM-6005a_1,7 of 100])
 FitResult(Sub[ADM-6005a_1,8 of 100])
 FitResult(Sub[ADM-6005a_1,9 of 100])
 FitResult(Sub[ADM-6005a_1,10 of 100])
 ⋮
 FitResult(Sub[ADM-6005a_1,92 of 100])
 FitResult(Sub[ADM-6005a_1,93 of 100])
 FitResult(Sub[ADM-6005a_1,94 of 100])
 FitResult(Sub[ADM-6005a_1,95 of 100])
 FitResult(Sub[ADM-6005a_1,96 of 100])
 FitResult(Sub[ADM-6005a_1,97 of 100])
 FitResult(Sub[ADM-6005a_1,98 of 100])
 FitResult(Sub[ADM-6005a_1,99 of 100])
 FitResult(Sub[ADM-6005a_1,100 of 100])

Summary statistics.

Repeat for the 15 measured spectra.

res= map(unks) do s
  fit_spectrum(s, frs)
end

15-element Vector{FilterFitResult{Float64}}:
 FitResult(ADM-6005a_1)
 FitResult(ADM-6005a_2)
 FitResult(ADM-6005a_3)
 FitResult(ADM-6005a_4)
 FitResult(ADM-6005a_5)
 FitResult(ADM-6005a_6)
 FitResult(ADM-6005a_7)
 FitResult(ADM-6005a_8)
 FitResult(ADM-6005a_9)
 FitResult(ADM-6005a_10)
 FitResult(ADM-6005a_11)
 FitResult(ADM-6005a_12)
 FitResult(ADM-6005a_13)
 FitResult(ADM-6005a_14)
 FitResult(ADM-6005a_15)

Summary statistics.

@btime fit_spectrum(unks[1], frs)

410.900 μs (1043 allocations: 611.44 KiB)
FitResult(ADM-6005a_1)