XAS standard deviation

The derivation of the standard deviation of the transmission, in XAS, is not consistent with the derivation of the mean transmission. The transmission is a weighted average with weights = Io, therefore the standard deviation should be the standard deviation of the weighted average.

See https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Weighted_sample_variance

We usually deal with large number of points, therefore the standard deviation is derived from the biased weighted sample variance.

The derivation of unbiased estimate of sample variance (reliablility weights) is also possibile but commented out, since it is more computationally demanding.

I checked in a few cases that both derivations give very similar results.

XAS.py

@@ -47,7 +47,6 @@ def absorption(T, Io, fluorescence=False):
     good = np.logical_and(np.isfinite(T), np.isfinite(Io))
     T = T[good]
     Io = Io[good]
-    
     # return type of the structured array
     fdtype = [('muA', 'f8'), ('sigmaA', 'f8'), ('weights', 'f8'),
         ('muT', 'f8'), ('sigmaT', 'f8'), ('muIo', 'f8'), ('sigmaIo', 'f8'),
@@ -64,21 +63,26 @@ def absorption(T, Io, fluorescence=False):
     weights = np.sum(Io)
 
     p = np.corrcoef(T, Io)[0,1]
-
-    if fluorescence:
-        muA = muT/muIo
-        sigmaA = np.abs(muA)*(np.sqrt((sigmaT/muT)**2 +
-            (sigmaIo/muIo)**2 - 2*p*sigmaIo*sigmaT/(muIo*muT)))
-    else:
-        muA = -np.log(muT/muIo)
-        # from error propagation for correlated data
-        sigmaA = (np.sqrt((sigmaT/muT)**2 + (sigmaIo/muIo)**2
-            - 2*p*sigmaIo*sigmaT/(muIo*muT)))
     
-        # direct calculation
-        #mask = (Io != 0)
-        #sigmaA = np.nanstd(-np.log(T[mask]/Io[mask]))
+    # weighted average of T/Io with Io as weights
+    muA = muT/muIo
+    
+    #Derivation of standard deviation
+    #1. using biased weighted sample variance:
+    #sigmaA = np.sqrt(np.average((T/Io - muA)**2, weights=Io))
 
+    #2. using unbiased weighted sample variance (reliablility weights):
+    V2 = np.sum(Io**2)
+    sigmaA = np.sqrt(np.sum(Io*(T/Io - muA)**2) / (weights - V2/weights))
+
+    #3. using error propagation for correlated data:
+    #sigmaA = np.abs(muA)*(np.sqrt((sigmaT/muT)**2 +
+    #    (sigmaIo/muIo)**2 - 2*p*sigmaIo*sigmaT/(muIo*muT)))
+    
+    if not fluorescence:
+        sigmaA = sigmaA/np.abs(muA)
+        muA = -np.log(muA)
+    
     return np.array([(muA, sigmaA, weights, muT, sigmaT, muIo, sigmaIo,
         p, counts)], dtype=fdtype)
 
@@ -97,7 +101,6 @@ def binning(x, data, func, bins=100, bin_length=None):
             bins: the bins edges
             res: a structured array of binned data
     """
-
     if bin_length is not None:
         bin_start = np.amin(x)
         bin_end = np.amax(x)
@@ -108,7 +111,7 @@ def binning(x, data, func, bins=100, bin_length=None):
         bins = np.linspace(bin_start, bin_end, bins)
     bin_centers = (bins[1:]+bins[:-1])/2
     nb_bins = len(bin_centers)
-
+    
     bin_idx = np.digitize(x, bins)
     dummy = func([])
     res = np.empty((nb_bins), dtype=dummy.dtype)
@@ -187,7 +190,7 @@ def xas(nrun, bins=None, Iokey='SCS_SA3', Itkey='MCP3apd', nrjkey='nrj', Iooffse
     dummy, nosample = binning(rundata['nrj'], rundata, whichIo, bins)
     muA = nosample['muA']
     sterrA = nosample['sigmaA']/np.sqrt(nosample['counts'])
-
+    
     bins_c = 0.5*(bins[1:] + bins[:-1])
     
     if plot: