added fitting errors to 1D scintillation parameter estimates

CHANGELOG.md

@@ -7,7 +7,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
-- `DynamicSpectrum` now copies colormaps to remove `MatplotlibDeprecationWarning`
+### Added
+
+- Fitting errors on the 1D scintillation parameters are now added. This changes the underlying `functionfit` output but this is not part of the public API.
+
+### Changed
+
+- `DynamicSpectrum` now copies colormaps with `copy.copy` to remove `MatplotlibDeprecationWarning`
 - Fixed scrunching and waterfall plot limits due to updates in numpy in the handling of NaNs. In waterfall plot, duplicate calls to set axis limits have been removed.
 
 ## [0.1.0] - 2021-05-03

pypulse/dynamicspectrum.py

@@ -237,20 +237,23 @@ def scintillation_parameters(self, plotbound=1.0, maxr=None, maxc=None,
             Taxis = (np.arange(-(NT-1), NT, dtype=np.float)*np.abs(dT))[1:-1] #???
 
 
-            pout, errs = ffit.gaussianfit(Taxis[NT//2:3*NT//2], self.acf[centerrind, NT//2:3*NT//2], baseline=True)
+            pout, perrs = ffit.gaussianfit(Taxis[NT//2:3*NT//2], self.acf[centerrind, NT//2:3*NT//2], baseline=True)
             ft = interpolate.interp1d(Taxis, ffit.funcgaussian(pout, Taxis, baseline=True)-(pout[3]+pout[0]/np.e))
             try:
                 delta_t_d = optimize.brentq(ft, 0, Taxis[-1])
+                err_t_d = perrs[2]
             except ValueError:
                 delta_t_d = np.nan
+                err_t_d = 0.0
 
-            pout, errs = ffit.gaussianfit(Faxis[NF//2:3*NF//2], self.acf[NF//2:3*NF//2, centercind], baseline=True)
+            pout, perrs = ffit.gaussianfit(Faxis[NF//2:3*NF//2], self.acf[NF//2:3*NF//2, centercind], baseline=True)
             fnu = interpolate.interp1d(Faxis, ffit.funcgaussian(pout, Faxis, baseline=True)-(pout[3]+pout[0]/2))
             try:
                 delta_nu_d = optimize.brentq(fnu, 0, Faxis[-1])
+                err_nu_d = perrs[2]
             except ValueError:
                 delta_nu_d = np.nan
-
+                err_t_d = 0.0
 
             #following Glenn's code and Cordes 1986 (Space Velocities...)
             # Errors from finite scintle effect:
@@ -265,8 +268,8 @@ def scintillation_parameters(self, plotbound=1.0, maxr=None, maxc=None,
             fse_nu_d = delta_nu_d/(2*np.log(2)*np.sqrt(N_d)) #log because of FWHM?
             fse_t_d = delta_t_d/(2*np.sqrt(N_d))
 
-            err_nu_d = fse_nu_d
-            err_t_d = fse_t_d #need to add in fitting errors
+            err_nu_d = np.sqrt(fse_nu_d**2 + err_nu_d**2)
+            err_t_d = np.sqrt(fse_t_d**2 + err_t_d**2)
 
             if show or savefig is not None:
                 fig = plt.figure()

pypulse/functionfit.py

@@ -21,10 +21,17 @@ def gaussianfit(x, y, baseline=False):
         p0 = [height, mu, sigma, 0.0]
     else:
         p0 = [height, mu, sigma]
-    p1, success = optimize.leastsq(errgaussian, p0[:], args=(x, y, baseline))
+    out = optimize.leastsq(errgaussian, p0[:], args=(x, y, baseline), full_output=1)
+    p1 = out[0] # output parameters
     p1[2] = np.abs(p1[2]) #enforces positive sigma
-    #Return values are the coefficients, the residuals
-    return p1, errgaussian(p1, x, y, baseline)
+
+    s_sq = (errgaussian(p1, x, y, baseline)**2).sum() / (len(x) - 4)
+    perrs = np.zeros_like(p1)
+    for i in range(len(perrs)):
+        perrs[i] = np.sqrt(out[1][i, i]*s_sq)
+
+    #Return values are the coefficients and the errors
+    return p1, perrs
 
 #area = np.sum(binwidths*hist)
 def funcsimpleDISSpdf(p, x, area=None):