From 12c2af50502bf75f3dcf6ea2b3e965985a206fc3 Mon Sep 17 00:00:00 2001 From: Michael Lam Date: Sat, 5 Jun 2021 20:34:37 -0400 Subject: [PATCH] added fitting errors to 1D scintillation parameter estimates --- CHANGELOG.md | 8 +++++++- pypulse/dynamicspectrum.py | 13 ++++++++----- pypulse/functionfit.py | 13 ++++++++++--- 3 files changed, 25 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c0a8d69..34a0882 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/pypulse/dynamicspectrum.py b/pypulse/dynamicspectrum.py index a9b6f44..6332f30 100644 --- a/pypulse/dynamicspectrum.py +++ b/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() diff --git a/pypulse/functionfit.py b/pypulse/functionfit.py index 2727b36..6a82387 100644 --- a/pypulse/functionfit.py +++ b/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):