Merge branch 'xas_error' into 'master'

XAS standard deviation

See merge request !80

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',
     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])
     bin_idx = np.digitize(rundata['nrj'], bins)
     nrj = np.array([ np.mean(rundata['nrj'][idx+1==bin_idx]) for idx in range(len(bins)-1)])