pyerrors.fits
1import gc 2from collections.abc import Sequence 3import warnings 4import numpy as np 5import autograd.numpy as anp 6import scipy.optimize 7import scipy.stats 8import matplotlib.pyplot as plt 9from matplotlib import gridspec 10from scipy.odr import ODR, Model, RealData 11from scipy.stats import chi2 12import iminuit 13from autograd import jacobian 14from autograd import elementwise_grad as egrad 15from .obs import Obs, derived_observable, covariance, cov_Obs 16 17 18class Fit_result(Sequence): 19 """Represents fit results. 20 21 Attributes 22 ---------- 23 fit_parameters : list 24 results for the individual fit parameters, 25 also accessible via indices. 26 """ 27 28 def __init__(self): 29 self.fit_parameters = None 30 31 def __getitem__(self, idx): 32 return self.fit_parameters[idx] 33 34 def __len__(self): 35 return len(self.fit_parameters) 36 37 def gamma_method(self): 38 """Apply the gamma method to all fit parameters""" 39 [o.gamma_method() for o in self.fit_parameters] 40 41 def __str__(self): 42 my_str = 'Goodness of fit:\n' 43 if hasattr(self, 'chisquare_by_dof'): 44 my_str += '\u03C7\u00b2/d.o.f. = ' + f'{self.chisquare_by_dof:2.6f}' + '\n' 45 elif hasattr(self, 'residual_variance'): 46 my_str += 'residual variance = ' + f'{self.residual_variance:2.6f}' + '\n' 47 if hasattr(self, 'chisquare_by_expected_chisquare'): 48 my_str += '\u03C7\u00b2/\u03C7\u00b2exp = ' + f'{self.chisquare_by_expected_chisquare:2.6f}' + '\n' 49 if hasattr(self, 'p_value'): 50 my_str += 'p-value = ' + f'{self.p_value:2.4f}' + '\n' 51 my_str += 'Fit parameters:\n' 52 for i_par, par in enumerate(self.fit_parameters): 53 my_str += str(i_par) + '\t' + ' ' * int(par >= 0) + str(par).rjust(int(par < 0.0)) + '\n' 54 return my_str 55 56 def __repr__(self): 57 m = max(map(len, list(self.__dict__.keys()))) + 1 58 return '\n'.join([key.rjust(m) + ': ' + repr(value) for key, value in sorted(self.__dict__.items())]) 59 60 61def least_squares(x, y, func, priors=None, silent=False, **kwargs): 62 r'''Performs a non-linear fit to y = func(x). 63 64 Parameters 65 ---------- 66 x : list 67 list of floats. 68 y : list 69 list of Obs. 70 func : object 71 fit function, has to be of the form 72 73 ```python 74 import autograd.numpy as anp 75 76 def func(a, x): 77 return a[0] + a[1] * x + a[2] * anp.sinh(x) 78 ``` 79 80 For multiple x values func can be of the form 81 82 ```python 83 def func(a, x): 84 (x1, x2) = x 85 return a[0] * x1 ** 2 + a[1] * x2 86 ``` 87 88 It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation 89 will not work. 90 priors : list, optional 91 priors has to be a list with an entry for every parameter in the fit. The entries can either be 92 Obs (e.g. results from a previous fit) or strings containing a value and an error formatted like 93 0.548(23), 500(40) or 0.5(0.4) 94 silent : bool, optional 95 If true all output to the console is omitted (default False). 96 initial_guess : list 97 can provide an initial guess for the input parameters. Relevant for 98 non-linear fits with many parameters. 99 method : str, optional 100 can be used to choose an alternative method for the minimization of chisquare. 101 The possible methods are the ones which can be used for scipy.optimize.minimize and 102 migrad of iminuit. If no method is specified, Levenberg-Marquard is used. 103 Reliable alternatives are migrad, Powell and Nelder-Mead. 104 correlated_fit : bool 105 If True, use the full inverse covariance matrix in the definition of the chisquare cost function. 106 For details about how the covariance matrix is estimated see `pyerrors.obs.covariance`. 107 In practice the correlation matrix is Cholesky decomposed and inverted (instead of the covariance matrix). 108 This procedure should be numerically more stable as the correlation matrix is typically better conditioned (Jacobi preconditioning). 109 At the moment this option only works for `prior==None` and when no `method` is given. 110 expected_chisquare : bool 111 If True estimates the expected chisquare which is 112 corrected by effects caused by correlated input data (default False). 113 resplot : bool 114 If True, a plot which displays fit, data and residuals is generated (default False). 115 qqplot : bool 116 If True, a quantile-quantile plot of the fit result is generated (default False). 117 ''' 118 if priors is not None: 119 return _prior_fit(x, y, func, priors, silent=silent, **kwargs) 120 else: 121 return _standard_fit(x, y, func, silent=silent, **kwargs) 122 123 124def total_least_squares(x, y, func, silent=False, **kwargs): 125 r'''Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters. 126 127 Parameters 128 ---------- 129 x : list 130 list of Obs, or a tuple of lists of Obs 131 y : list 132 list of Obs. The dvalues of the Obs are used as x- and yerror for the fit. 133 func : object 134 func has to be of the form 135 136 ```python 137 import autograd.numpy as anp 138 139 def func(a, x): 140 return a[0] + a[1] * x + a[2] * anp.sinh(x) 141 ``` 142 143 For multiple x values func can be of the form 144 145 ```python 146 def func(a, x): 147 (x1, x2) = x 148 return a[0] * x1 ** 2 + a[1] * x2 149 ``` 150 151 It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation 152 will not work. 153 silent : bool, optional 154 If true all output to the console is omitted (default False). 155 initial_guess : list 156 can provide an initial guess for the input parameters. Relevant for non-linear 157 fits with many parameters. 158 expected_chisquare : bool 159 If true prints the expected chisquare which is 160 corrected by effects caused by correlated input data. 161 This can take a while as the full correlation matrix 162 has to be calculated (default False). 163 164 Notes 165 ----- 166 Based on the orthogonal distance regression module of scipy 167 ''' 168 169 output = Fit_result() 170 171 output.fit_function = func 172 173 x = np.array(x) 174 175 x_shape = x.shape 176 177 if not callable(func): 178 raise TypeError('func has to be a function.') 179 180 for i in range(25): 181 try: 182 func(np.arange(i), x.T[0]) 183 except Exception: 184 pass 185 else: 186 break 187 188 n_parms = i 189 if not silent: 190 print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1)) 191 192 x_f = np.vectorize(lambda o: o.value)(x) 193 dx_f = np.vectorize(lambda o: o.dvalue)(x) 194 y_f = np.array([o.value for o in y]) 195 dy_f = np.array([o.dvalue for o in y]) 196 197 if np.any(np.asarray(dx_f) <= 0.0): 198 raise Exception('No x errors available, run the gamma method first.') 199 200 if np.any(np.asarray(dy_f) <= 0.0): 201 raise Exception('No y errors available, run the gamma method first.') 202 203 if 'initial_guess' in kwargs: 204 x0 = kwargs.get('initial_guess') 205 if len(x0) != n_parms: 206 raise Exception('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms)) 207 else: 208 x0 = [1] * n_parms 209 210 data = RealData(x_f, y_f, sx=dx_f, sy=dy_f) 211 model = Model(func) 212 odr = ODR(data, model, x0, partol=np.finfo(np.float64).eps) 213 odr.set_job(fit_type=0, deriv=1) 214 out = odr.run() 215 216 output.residual_variance = out.res_var 217 218 output.method = 'ODR' 219 220 output.message = out.stopreason 221 222 output.xplus = out.xplus 223 224 if not silent: 225 print('Method: ODR') 226 print(*out.stopreason) 227 print('Residual variance:', output.residual_variance) 228 229 if out.info > 3: 230 raise Exception('The minimization procedure did not converge.') 231 232 m = x_f.size 233 234 def odr_chisquare(p): 235 model = func(p[:n_parms], p[n_parms:].reshape(x_shape)) 236 chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((x_f - p[n_parms:].reshape(x_shape)) / dx_f) ** 2) 237 return chisq 238 239 if kwargs.get('expected_chisquare') is True: 240 W = np.diag(1 / np.asarray(np.concatenate((dy_f.ravel(), dx_f.ravel())))) 241 242 if kwargs.get('covariance') is not None: 243 cov = kwargs.get('covariance') 244 else: 245 cov = covariance(np.concatenate((y, x.ravel()))) 246 247 number_of_x_parameters = int(m / x_f.shape[-1]) 248 249 old_jac = jacobian(func)(out.beta, out.xplus) 250 fused_row1 = np.concatenate((old_jac, np.concatenate((number_of_x_parameters * [np.zeros(old_jac.shape)]), axis=0))) 251 fused_row2 = np.concatenate((jacobian(lambda x, y: func(y, x))(out.xplus, out.beta).reshape(x_f.shape[-1], x_f.shape[-1] * number_of_x_parameters), np.identity(number_of_x_parameters * old_jac.shape[0]))) 252 new_jac = np.concatenate((fused_row1, fused_row2), axis=1) 253 254 A = W @ new_jac 255 P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T 256 expected_chisquare = np.trace((np.identity(P_phi.shape[0]) - P_phi) @ W @ cov @ W) 257 if expected_chisquare <= 0.0: 258 warnings.warn("Negative expected_chisquare.", RuntimeWarning) 259 expected_chisquare = np.abs(expected_chisquare) 260 output.chisquare_by_expected_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel()))) / expected_chisquare 261 if not silent: 262 print('chisquare/expected_chisquare:', 263 output.chisquare_by_expected_chisquare) 264 265 fitp = out.beta 266 try: 267 hess = jacobian(jacobian(odr_chisquare))(np.concatenate((fitp, out.xplus.ravel()))) 268 except TypeError: 269 raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None 270 271 def odr_chisquare_compact_x(d): 272 model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape)) 273 chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((d[n_parms + m:].reshape(x_shape) - d[n_parms:n_parms + m].reshape(x_shape)) / dx_f) ** 2) 274 return chisq 275 276 jac_jac_x = jacobian(jacobian(odr_chisquare_compact_x))(np.concatenate((fitp, out.xplus.ravel(), x_f.ravel()))) 277 278 # Compute hess^{-1} @ jac_jac_x[:n_parms + m, n_parms + m:] using LAPACK dgesv 279 try: 280 deriv_x = -scipy.linalg.solve(hess, jac_jac_x[:n_parms + m, n_parms + m:]) 281 except np.linalg.LinAlgError: 282 raise Exception("Cannot invert hessian matrix.") 283 284 def odr_chisquare_compact_y(d): 285 model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape)) 286 chisq = anp.sum(((d[n_parms + m:] - model) / dy_f) ** 2) + anp.sum(((x_f - d[n_parms:n_parms + m].reshape(x_shape)) / dx_f) ** 2) 287 return chisq 288 289 jac_jac_y = jacobian(jacobian(odr_chisquare_compact_y))(np.concatenate((fitp, out.xplus.ravel(), y_f))) 290 291 # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv 292 try: 293 deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms + m, n_parms + m:]) 294 except np.linalg.LinAlgError: 295 raise Exception("Cannot invert hessian matrix.") 296 297 result = [] 298 for i in range(n_parms): 299 result.append(derived_observable(lambda my_var, **kwargs: (my_var[0] + np.finfo(np.float64).eps) / (x.ravel()[0].value + np.finfo(np.float64).eps) * out.beta[i], list(x.ravel()) + list(y), man_grad=list(deriv_x[i]) + list(deriv_y[i]))) 300 301 output.fit_parameters = result 302 303 output.odr_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel()))) 304 output.dof = x.shape[-1] - n_parms 305 output.p_value = 1 - chi2.cdf(output.odr_chisquare, output.dof) 306 307 return output 308 309 310def _prior_fit(x, y, func, priors, silent=False, **kwargs): 311 output = Fit_result() 312 313 output.fit_function = func 314 315 x = np.asarray(x) 316 317 if not callable(func): 318 raise TypeError('func has to be a function.') 319 320 for i in range(100): 321 try: 322 func(np.arange(i), 0) 323 except Exception: 324 pass 325 else: 326 break 327 328 n_parms = i 329 330 if n_parms != len(priors): 331 raise Exception('Priors does not have the correct length.') 332 333 def extract_val_and_dval(string): 334 split_string = string.split('(') 335 if '.' in split_string[0] and '.' not in split_string[1][:-1]: 336 factor = 10 ** -len(split_string[0].partition('.')[2]) 337 else: 338 factor = 1 339 return float(split_string[0]), float(split_string[1][:-1]) * factor 340 341 loc_priors = [] 342 for i_n, i_prior in enumerate(priors): 343 if isinstance(i_prior, Obs): 344 loc_priors.append(i_prior) 345 else: 346 loc_val, loc_dval = extract_val_and_dval(i_prior) 347 loc_priors.append(cov_Obs(loc_val, loc_dval ** 2, '#prior' + str(i_n) + f"_{np.random.randint(2147483647):010d}")) 348 349 output.priors = loc_priors 350 351 if not silent: 352 print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1)) 353 354 y_f = [o.value for o in y] 355 dy_f = [o.dvalue for o in y] 356 357 if np.any(np.asarray(dy_f) <= 0.0): 358 raise Exception('No y errors available, run the gamma method first.') 359 360 p_f = [o.value for o in loc_priors] 361 dp_f = [o.dvalue for o in loc_priors] 362 363 if np.any(np.asarray(dp_f) <= 0.0): 364 raise Exception('No prior errors available, run the gamma method first.') 365 366 if 'initial_guess' in kwargs: 367 x0 = kwargs.get('initial_guess') 368 if len(x0) != n_parms: 369 raise Exception('Initial guess does not have the correct length.') 370 else: 371 x0 = p_f 372 373 def chisqfunc(p): 374 model = func(p, x) 375 chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((p_f - p) / dp_f) ** 2) 376 return chisq 377 378 if not silent: 379 print('Method: migrad') 380 381 m = iminuit.Minuit(chisqfunc, x0) 382 m.errordef = 1 383 m.print_level = 0 384 if 'tol' in kwargs: 385 m.tol = kwargs.get('tol') 386 else: 387 m.tol = 1e-4 388 m.migrad() 389 params = np.asarray(m.values) 390 391 output.chisquare_by_dof = m.fval / len(x) 392 393 output.method = 'migrad' 394 395 if not silent: 396 print('chisquare/d.o.f.:', output.chisquare_by_dof) 397 398 if not m.fmin.is_valid: 399 raise Exception('The minimization procedure did not converge.') 400 401 hess_inv = np.linalg.pinv(jacobian(jacobian(chisqfunc))(params)) 402 403 def chisqfunc_compact(d): 404 model = func(d[:n_parms], x) 405 chisq = anp.sum(((d[n_parms: n_parms + len(x)] - model) / dy_f) ** 2) + anp.sum(((d[n_parms + len(x):] - d[:n_parms]) / dp_f) ** 2) 406 return chisq 407 408 jac_jac = jacobian(jacobian(chisqfunc_compact))(np.concatenate((params, y_f, p_f))) 409 410 deriv = -hess_inv @ jac_jac[:n_parms, n_parms:] 411 412 result = [] 413 for i in range(n_parms): 414 result.append(derived_observable(lambda x, **kwargs: (x[0] + np.finfo(np.float64).eps) / (y[0].value + np.finfo(np.float64).eps) * params[i], list(y) + list(loc_priors), man_grad=list(deriv[i]))) 415 416 output.fit_parameters = result 417 output.chisquare = chisqfunc(np.asarray(params)) 418 419 if kwargs.get('resplot') is True: 420 residual_plot(x, y, func, result) 421 422 if kwargs.get('qqplot') is True: 423 qqplot(x, y, func, result) 424 425 return output 426 427 428def _standard_fit(x, y, func, silent=False, **kwargs): 429 430 output = Fit_result() 431 432 output.fit_function = func 433 434 x = np.asarray(x) 435 436 if x.shape[-1] != len(y): 437 raise Exception('x and y input have to have the same length') 438 439 if len(x.shape) > 2: 440 raise Exception('Unknown format for x values') 441 442 if not callable(func): 443 raise TypeError('func has to be a function.') 444 445 for i in range(25): 446 try: 447 func(np.arange(i), x.T[0]) 448 except Exception: 449 pass 450 else: 451 break 452 453 n_parms = i 454 455 if not silent: 456 print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1)) 457 458 y_f = [o.value for o in y] 459 dy_f = [o.dvalue for o in y] 460 461 if np.any(np.asarray(dy_f) <= 0.0): 462 raise Exception('No y errors available, run the gamma method first.') 463 464 if 'initial_guess' in kwargs: 465 x0 = kwargs.get('initial_guess') 466 if len(x0) != n_parms: 467 raise Exception('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms)) 468 else: 469 x0 = [0.1] * n_parms 470 471 if kwargs.get('correlated_fit') is True: 472 corr = covariance(y, correlation=True, **kwargs) 473 covdiag = np.diag(1 / np.asarray(dy_f)) 474 condn = np.linalg.cond(corr) 475 if condn > 0.1 / np.finfo(float).eps: 476 raise Exception(f"Cannot invert correlation matrix as its condition number exceeds machine precision ({condn:1.2e})") 477 if condn > 1 / np.sqrt(np.finfo(float).eps): 478 warnings.warn("Correlation matrix may be ill-conditioned, condition number: {%1.2e}" % (condn), RuntimeWarning) 479 chol = np.linalg.cholesky(corr) 480 chol_inv = scipy.linalg.solve(chol, covdiag) 481 482 def chisqfunc_corr(p): 483 model = func(p, x) 484 chisq = anp.sum(anp.dot(chol_inv, (y_f - model)) ** 2) 485 return chisq 486 487 def chisqfunc(p): 488 model = func(p, x) 489 chisq = anp.sum(((y_f - model) / dy_f) ** 2) 490 return chisq 491 492 output.method = kwargs.get('method', 'Levenberg-Marquardt') 493 if not silent: 494 print('Method:', output.method) 495 496 if output.method != 'Levenberg-Marquardt': 497 if output.method == 'migrad': 498 fit_result = iminuit.minimize(chisqfunc, x0, tol=1e-4) # Stopping criterion 0.002 * tol * errordef 499 if kwargs.get('correlated_fit') is True: 500 fit_result = iminuit.minimize(chisqfunc_corr, fit_result.x, tol=1e-4) # Stopping criterion 0.002 * tol * errordef 501 output.iterations = fit_result.nfev 502 else: 503 fit_result = scipy.optimize.minimize(chisqfunc, x0, method=kwargs.get('method'), tol=1e-12) 504 if kwargs.get('correlated_fit') is True: 505 fit_result = scipy.optimize.minimize(chisqfunc_corr, fit_result.x, method=kwargs.get('method'), tol=1e-12) 506 output.iterations = fit_result.nit 507 508 chisquare = fit_result.fun 509 510 else: 511 if kwargs.get('correlated_fit') is True: 512 def chisqfunc_residuals_corr(p): 513 model = func(p, x) 514 chisq = anp.dot(chol_inv, (y_f - model)) 515 return chisq 516 517 def chisqfunc_residuals(p): 518 model = func(p, x) 519 chisq = ((y_f - model) / dy_f) 520 return chisq 521 522 fit_result = scipy.optimize.least_squares(chisqfunc_residuals, x0, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15) 523 if kwargs.get('correlated_fit') is True: 524 fit_result = scipy.optimize.least_squares(chisqfunc_residuals_corr, fit_result.x, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15) 525 526 chisquare = np.sum(fit_result.fun ** 2) 527 if kwargs.get('correlated_fit') is True: 528 assert np.isclose(chisquare, chisqfunc_corr(fit_result.x), atol=1e-14) 529 else: 530 assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14) 531 532 output.iterations = fit_result.nfev 533 534 if not fit_result.success: 535 raise Exception('The minimization procedure did not converge.') 536 537 if x.shape[-1] - n_parms > 0: 538 output.chisquare_by_dof = chisquare / (x.shape[-1] - n_parms) 539 else: 540 output.chisquare_by_dof = float('nan') 541 542 output.message = fit_result.message 543 if not silent: 544 print(fit_result.message) 545 print('chisquare/d.o.f.:', output.chisquare_by_dof) 546 547 if kwargs.get('expected_chisquare') is True: 548 if kwargs.get('correlated_fit') is not True: 549 W = np.diag(1 / np.asarray(dy_f)) 550 cov = covariance(y) 551 A = W @ jacobian(func)(fit_result.x, x) 552 P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T 553 expected_chisquare = np.trace((np.identity(x.shape[-1]) - P_phi) @ W @ cov @ W) 554 output.chisquare_by_expected_chisquare = chisquare / expected_chisquare 555 if not silent: 556 print('chisquare/expected_chisquare:', 557 output.chisquare_by_expected_chisquare) 558 559 fitp = fit_result.x 560 try: 561 hess = jacobian(jacobian(chisqfunc))(fitp) 562 except TypeError: 563 raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None 564 565 if kwargs.get('correlated_fit') is True: 566 def chisqfunc_compact(d): 567 model = func(d[:n_parms], x) 568 chisq = anp.sum(anp.dot(chol_inv, (d[n_parms:] - model)) ** 2) 569 return chisq 570 571 else: 572 def chisqfunc_compact(d): 573 model = func(d[:n_parms], x) 574 chisq = anp.sum(((d[n_parms:] - model) / dy_f) ** 2) 575 return chisq 576 577 jac_jac = jacobian(jacobian(chisqfunc_compact))(np.concatenate((fitp, y_f))) 578 579 # Compute hess^{-1} @ jac_jac[:n_parms, n_parms:] using LAPACK dgesv 580 try: 581 deriv = -scipy.linalg.solve(hess, jac_jac[:n_parms, n_parms:]) 582 except np.linalg.LinAlgError: 583 raise Exception("Cannot invert hessian matrix.") 584 585 result = [] 586 for i in range(n_parms): 587 result.append(derived_observable(lambda x, **kwargs: (x[0] + np.finfo(np.float64).eps) / (y[0].value + np.finfo(np.float64).eps) * fit_result.x[i], list(y), man_grad=list(deriv[i]))) 588 589 output.fit_parameters = result 590 591 output.chisquare = chisquare 592 output.dof = x.shape[-1] - n_parms 593 output.p_value = 1 - chi2.cdf(output.chisquare, output.dof) 594 595 if kwargs.get('resplot') is True: 596 residual_plot(x, y, func, result) 597 598 if kwargs.get('qqplot') is True: 599 qqplot(x, y, func, result) 600 601 return output 602 603 604def fit_lin(x, y, **kwargs): 605 """Performs a linear fit to y = n + m * x and returns two Obs n, m. 606 607 Parameters 608 ---------- 609 x : list 610 Can either be a list of floats in which case no xerror is assumed, or 611 a list of Obs, where the dvalues of the Obs are used as xerror for the fit. 612 y : list 613 List of Obs, the dvalues of the Obs are used as yerror for the fit. 614 """ 615 616 def f(a, x): 617 y = a[0] + a[1] * x 618 return y 619 620 if all(isinstance(n, Obs) for n in x): 621 out = total_least_squares(x, y, f, **kwargs) 622 return out.fit_parameters 623 elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray): 624 out = least_squares(x, y, f, **kwargs) 625 return out.fit_parameters 626 else: 627 raise Exception('Unsupported types for x') 628 629 630def qqplot(x, o_y, func, p): 631 """Generates a quantile-quantile plot of the fit result which can be used to 632 check if the residuals of the fit are gaussian distributed. 633 """ 634 635 residuals = [] 636 for i_x, i_y in zip(x, o_y): 637 residuals.append((i_y - func(p, i_x)) / i_y.dvalue) 638 residuals = sorted(residuals) 639 my_y = [o.value for o in residuals] 640 probplot = scipy.stats.probplot(my_y) 641 my_x = probplot[0][0] 642 plt.figure(figsize=(8, 8 / 1.618)) 643 plt.errorbar(my_x, my_y, fmt='o') 644 fit_start = my_x[0] 645 fit_stop = my_x[-1] 646 samples = np.arange(fit_start, fit_stop, 0.01) 647 plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution') 648 plt.plot(samples, probplot[1][0] * samples + probplot[1][1], zorder=10, label='Least squares fit, r=' + str(np.around(probplot[1][2], 3)), marker='', ls='-') 649 650 plt.xlabel('Theoretical quantiles') 651 plt.ylabel('Ordered Values') 652 plt.legend() 653 plt.draw() 654 655 656def residual_plot(x, y, func, fit_res): 657 """ Generates a plot which compares the fit to the data and displays the corresponding residuals""" 658 sorted_x = sorted(x) 659 xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0]) 660 xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2]) 661 x_samples = np.arange(xstart, xstop + 0.01, 0.01) 662 663 plt.figure(figsize=(8, 8 / 1.618)) 664 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0) 665 ax0 = plt.subplot(gs[0]) 666 ax0.errorbar(x, [o.value for o in y], yerr=[o.dvalue for o in y], ls='none', fmt='o', capsize=3, markersize=5, label='Data') 667 ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0) 668 ax0.set_xticklabels([]) 669 ax0.set_xlim([xstart, xstop]) 670 ax0.set_xticklabels([]) 671 ax0.legend() 672 673 residuals = (np.asarray([o.value for o in y]) - func([o.value for o in fit_res], x)) / np.asarray([o.dvalue for o in y]) 674 ax1 = plt.subplot(gs[1]) 675 ax1.plot(x, residuals, 'ko', ls='none', markersize=5) 676 ax1.tick_params(direction='out') 677 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True) 678 ax1.axhline(y=0.0, ls='--', color='k', marker=" ") 679 ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k') 680 ax1.set_xlim([xstart, xstop]) 681 ax1.set_ylabel('Residuals') 682 plt.subplots_adjust(wspace=None, hspace=None) 683 plt.draw() 684 685 686def error_band(x, func, beta): 687 """Returns the error band for an array of sample values x, for given fit function func with optimized parameters beta.""" 688 cov = covariance(beta) 689 if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps): 690 warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning) 691 692 deriv = [] 693 for i, item in enumerate(x): 694 deriv.append(np.array(egrad(func)([o.value for o in beta], item))) 695 696 err = [] 697 for i, item in enumerate(x): 698 err.append(np.sqrt(deriv[i] @ cov @ deriv[i])) 699 err = np.array(err) 700 701 return err 702 703 704def ks_test(objects=None): 705 """Performs a Kolmogorov–Smirnov test for the p-values of all fit object. 706 707 Parameters 708 ---------- 709 objects : list 710 List of fit results to include in the analysis (optional). 711 """ 712 713 if objects is None: 714 obs_list = [] 715 for obj in gc.get_objects(): 716 if isinstance(obj, Fit_result): 717 obs_list.append(obj) 718 else: 719 obs_list = objects 720 721 p_values = [o.p_value for o in obs_list] 722 723 bins = len(p_values) 724 x = np.arange(0, 1.001, 0.001) 725 plt.plot(x, x, 'k', zorder=1) 726 plt.xlim(0, 1) 727 plt.ylim(0, 1) 728 plt.xlabel('p-value') 729 plt.ylabel('Cumulative probability') 730 plt.title(str(bins) + ' p-values') 731 732 n = np.arange(1, bins + 1) / np.float64(bins) 733 Xs = np.sort(p_values) 734 plt.step(Xs, n) 735 diffs = n - Xs 736 loc_max_diff = np.argmax(np.abs(diffs)) 737 loc = Xs[loc_max_diff] 738 plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0)) 739 plt.draw() 740 741 print(scipy.stats.kstest(p_values, 'uniform'))
19class Fit_result(Sequence): 20 """Represents fit results. 21 22 Attributes 23 ---------- 24 fit_parameters : list 25 results for the individual fit parameters, 26 also accessible via indices. 27 """ 28 29 def __init__(self): 30 self.fit_parameters = None 31 32 def __getitem__(self, idx): 33 return self.fit_parameters[idx] 34 35 def __len__(self): 36 return len(self.fit_parameters) 37 38 def gamma_method(self): 39 """Apply the gamma method to all fit parameters""" 40 [o.gamma_method() for o in self.fit_parameters] 41 42 def __str__(self): 43 my_str = 'Goodness of fit:\n' 44 if hasattr(self, 'chisquare_by_dof'): 45 my_str += '\u03C7\u00b2/d.o.f. = ' + f'{self.chisquare_by_dof:2.6f}' + '\n' 46 elif hasattr(self, 'residual_variance'): 47 my_str += 'residual variance = ' + f'{self.residual_variance:2.6f}' + '\n' 48 if hasattr(self, 'chisquare_by_expected_chisquare'): 49 my_str += '\u03C7\u00b2/\u03C7\u00b2exp = ' + f'{self.chisquare_by_expected_chisquare:2.6f}' + '\n' 50 if hasattr(self, 'p_value'): 51 my_str += 'p-value = ' + f'{self.p_value:2.4f}' + '\n' 52 my_str += 'Fit parameters:\n' 53 for i_par, par in enumerate(self.fit_parameters): 54 my_str += str(i_par) + '\t' + ' ' * int(par >= 0) + str(par).rjust(int(par < 0.0)) + '\n' 55 return my_str 56 57 def __repr__(self): 58 m = max(map(len, list(self.__dict__.keys()))) + 1 59 return '\n'.join([key.rjust(m) + ': ' + repr(value) for key, value in sorted(self.__dict__.items())])
Represents fit results.
Attributes
- fit_parameters (list): results for the individual fit parameters, also accessible via indices.
38 def gamma_method(self): 39 """Apply the gamma method to all fit parameters""" 40 [o.gamma_method() for o in self.fit_parameters]
Apply the gamma method to all fit parameters
Inherited Members
- collections.abc.Sequence
- index
- count
62def least_squares(x, y, func, priors=None, silent=False, **kwargs): 63 r'''Performs a non-linear fit to y = func(x). 64 65 Parameters 66 ---------- 67 x : list 68 list of floats. 69 y : list 70 list of Obs. 71 func : object 72 fit function, has to be of the form 73 74 ```python 75 import autograd.numpy as anp 76 77 def func(a, x): 78 return a[0] + a[1] * x + a[2] * anp.sinh(x) 79 ``` 80 81 For multiple x values func can be of the form 82 83 ```python 84 def func(a, x): 85 (x1, x2) = x 86 return a[0] * x1 ** 2 + a[1] * x2 87 ``` 88 89 It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation 90 will not work. 91 priors : list, optional 92 priors has to be a list with an entry for every parameter in the fit. The entries can either be 93 Obs (e.g. results from a previous fit) or strings containing a value and an error formatted like 94 0.548(23), 500(40) or 0.5(0.4) 95 silent : bool, optional 96 If true all output to the console is omitted (default False). 97 initial_guess : list 98 can provide an initial guess for the input parameters. Relevant for 99 non-linear fits with many parameters. 100 method : str, optional 101 can be used to choose an alternative method for the minimization of chisquare. 102 The possible methods are the ones which can be used for scipy.optimize.minimize and 103 migrad of iminuit. If no method is specified, Levenberg-Marquard is used. 104 Reliable alternatives are migrad, Powell and Nelder-Mead. 105 correlated_fit : bool 106 If True, use the full inverse covariance matrix in the definition of the chisquare cost function. 107 For details about how the covariance matrix is estimated see `pyerrors.obs.covariance`. 108 In practice the correlation matrix is Cholesky decomposed and inverted (instead of the covariance matrix). 109 This procedure should be numerically more stable as the correlation matrix is typically better conditioned (Jacobi preconditioning). 110 At the moment this option only works for `prior==None` and when no `method` is given. 111 expected_chisquare : bool 112 If True estimates the expected chisquare which is 113 corrected by effects caused by correlated input data (default False). 114 resplot : bool 115 If True, a plot which displays fit, data and residuals is generated (default False). 116 qqplot : bool 117 If True, a quantile-quantile plot of the fit result is generated (default False). 118 ''' 119 if priors is not None: 120 return _prior_fit(x, y, func, priors, silent=silent, **kwargs) 121 else: 122 return _standard_fit(x, y, func, silent=silent, **kwargs)
Performs a non-linear fit to y = func(x).
Parameters
- x (list): list of floats.
- y (list): list of Obs.
func (object): fit function, has to be of the form
import autograd.numpy as anp def func(a, x): return a[0] + a[1] * x + a[2] * anp.sinh(x)
For multiple x values func can be of the form
def func(a, x): (x1, x2) = x return a[0] * x1 ** 2 + a[1] * x2
It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation will not work.
- priors (list, optional): priors has to be a list with an entry for every parameter in the fit. The entries can either be Obs (e.g. results from a previous fit) or strings containing a value and an error formatted like 0.548(23), 500(40) or 0.5(0.4)
- silent (bool, optional): If true all output to the console is omitted (default False).
- initial_guess (list): can provide an initial guess for the input parameters. Relevant for non-linear fits with many parameters.
- method (str, optional): can be used to choose an alternative method for the minimization of chisquare. The possible methods are the ones which can be used for scipy.optimize.minimize and migrad of iminuit. If no method is specified, Levenberg-Marquard is used. Reliable alternatives are migrad, Powell and Nelder-Mead.
- correlated_fit (bool):
If True, use the full inverse covariance matrix in the definition of the chisquare cost function.
For details about how the covariance matrix is estimated see
pyerrors.obs.covariance
. In practice the correlation matrix is Cholesky decomposed and inverted (instead of the covariance matrix). This procedure should be numerically more stable as the correlation matrix is typically better conditioned (Jacobi preconditioning). At the moment this option only works forprior==None
and when nomethod
is given. - expected_chisquare (bool): If True estimates the expected chisquare which is corrected by effects caused by correlated input data (default False).
- resplot (bool): If True, a plot which displays fit, data and residuals is generated (default False).
- qqplot (bool): If True, a quantile-quantile plot of the fit result is generated (default False).
125def total_least_squares(x, y, func, silent=False, **kwargs): 126 r'''Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters. 127 128 Parameters 129 ---------- 130 x : list 131 list of Obs, or a tuple of lists of Obs 132 y : list 133 list of Obs. The dvalues of the Obs are used as x- and yerror for the fit. 134 func : object 135 func has to be of the form 136 137 ```python 138 import autograd.numpy as anp 139 140 def func(a, x): 141 return a[0] + a[1] * x + a[2] * anp.sinh(x) 142 ``` 143 144 For multiple x values func can be of the form 145 146 ```python 147 def func(a, x): 148 (x1, x2) = x 149 return a[0] * x1 ** 2 + a[1] * x2 150 ``` 151 152 It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation 153 will not work. 154 silent : bool, optional 155 If true all output to the console is omitted (default False). 156 initial_guess : list 157 can provide an initial guess for the input parameters. Relevant for non-linear 158 fits with many parameters. 159 expected_chisquare : bool 160 If true prints the expected chisquare which is 161 corrected by effects caused by correlated input data. 162 This can take a while as the full correlation matrix 163 has to be calculated (default False). 164 165 Notes 166 ----- 167 Based on the orthogonal distance regression module of scipy 168 ''' 169 170 output = Fit_result() 171 172 output.fit_function = func 173 174 x = np.array(x) 175 176 x_shape = x.shape 177 178 if not callable(func): 179 raise TypeError('func has to be a function.') 180 181 for i in range(25): 182 try: 183 func(np.arange(i), x.T[0]) 184 except Exception: 185 pass 186 else: 187 break 188 189 n_parms = i 190 if not silent: 191 print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1)) 192 193 x_f = np.vectorize(lambda o: o.value)(x) 194 dx_f = np.vectorize(lambda o: o.dvalue)(x) 195 y_f = np.array([o.value for o in y]) 196 dy_f = np.array([o.dvalue for o in y]) 197 198 if np.any(np.asarray(dx_f) <= 0.0): 199 raise Exception('No x errors available, run the gamma method first.') 200 201 if np.any(np.asarray(dy_f) <= 0.0): 202 raise Exception('No y errors available, run the gamma method first.') 203 204 if 'initial_guess' in kwargs: 205 x0 = kwargs.get('initial_guess') 206 if len(x0) != n_parms: 207 raise Exception('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms)) 208 else: 209 x0 = [1] * n_parms 210 211 data = RealData(x_f, y_f, sx=dx_f, sy=dy_f) 212 model = Model(func) 213 odr = ODR(data, model, x0, partol=np.finfo(np.float64).eps) 214 odr.set_job(fit_type=0, deriv=1) 215 out = odr.run() 216 217 output.residual_variance = out.res_var 218 219 output.method = 'ODR' 220 221 output.message = out.stopreason 222 223 output.xplus = out.xplus 224 225 if not silent: 226 print('Method: ODR') 227 print(*out.stopreason) 228 print('Residual variance:', output.residual_variance) 229 230 if out.info > 3: 231 raise Exception('The minimization procedure did not converge.') 232 233 m = x_f.size 234 235 def odr_chisquare(p): 236 model = func(p[:n_parms], p[n_parms:].reshape(x_shape)) 237 chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((x_f - p[n_parms:].reshape(x_shape)) / dx_f) ** 2) 238 return chisq 239 240 if kwargs.get('expected_chisquare') is True: 241 W = np.diag(1 / np.asarray(np.concatenate((dy_f.ravel(), dx_f.ravel())))) 242 243 if kwargs.get('covariance') is not None: 244 cov = kwargs.get('covariance') 245 else: 246 cov = covariance(np.concatenate((y, x.ravel()))) 247 248 number_of_x_parameters = int(m / x_f.shape[-1]) 249 250 old_jac = jacobian(func)(out.beta, out.xplus) 251 fused_row1 = np.concatenate((old_jac, np.concatenate((number_of_x_parameters * [np.zeros(old_jac.shape)]), axis=0))) 252 fused_row2 = np.concatenate((jacobian(lambda x, y: func(y, x))(out.xplus, out.beta).reshape(x_f.shape[-1], x_f.shape[-1] * number_of_x_parameters), np.identity(number_of_x_parameters * old_jac.shape[0]))) 253 new_jac = np.concatenate((fused_row1, fused_row2), axis=1) 254 255 A = W @ new_jac 256 P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T 257 expected_chisquare = np.trace((np.identity(P_phi.shape[0]) - P_phi) @ W @ cov @ W) 258 if expected_chisquare <= 0.0: 259 warnings.warn("Negative expected_chisquare.", RuntimeWarning) 260 expected_chisquare = np.abs(expected_chisquare) 261 output.chisquare_by_expected_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel()))) / expected_chisquare 262 if not silent: 263 print('chisquare/expected_chisquare:', 264 output.chisquare_by_expected_chisquare) 265 266 fitp = out.beta 267 try: 268 hess = jacobian(jacobian(odr_chisquare))(np.concatenate((fitp, out.xplus.ravel()))) 269 except TypeError: 270 raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None 271 272 def odr_chisquare_compact_x(d): 273 model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape)) 274 chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((d[n_parms + m:].reshape(x_shape) - d[n_parms:n_parms + m].reshape(x_shape)) / dx_f) ** 2) 275 return chisq 276 277 jac_jac_x = jacobian(jacobian(odr_chisquare_compact_x))(np.concatenate((fitp, out.xplus.ravel(), x_f.ravel()))) 278 279 # Compute hess^{-1} @ jac_jac_x[:n_parms + m, n_parms + m:] using LAPACK dgesv 280 try: 281 deriv_x = -scipy.linalg.solve(hess, jac_jac_x[:n_parms + m, n_parms + m:]) 282 except np.linalg.LinAlgError: 283 raise Exception("Cannot invert hessian matrix.") 284 285 def odr_chisquare_compact_y(d): 286 model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape)) 287 chisq = anp.sum(((d[n_parms + m:] - model) / dy_f) ** 2) + anp.sum(((x_f - d[n_parms:n_parms + m].reshape(x_shape)) / dx_f) ** 2) 288 return chisq 289 290 jac_jac_y = jacobian(jacobian(odr_chisquare_compact_y))(np.concatenate((fitp, out.xplus.ravel(), y_f))) 291 292 # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv 293 try: 294 deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms + m, n_parms + m:]) 295 except np.linalg.LinAlgError: 296 raise Exception("Cannot invert hessian matrix.") 297 298 result = [] 299 for i in range(n_parms): 300 result.append(derived_observable(lambda my_var, **kwargs: (my_var[0] + np.finfo(np.float64).eps) / (x.ravel()[0].value + np.finfo(np.float64).eps) * out.beta[i], list(x.ravel()) + list(y), man_grad=list(deriv_x[i]) + list(deriv_y[i]))) 301 302 output.fit_parameters = result 303 304 output.odr_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel()))) 305 output.dof = x.shape[-1] - n_parms 306 output.p_value = 1 - chi2.cdf(output.odr_chisquare, output.dof) 307 308 return output
Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters.
Parameters
- x (list): list of Obs, or a tuple of lists of Obs
- y (list): list of Obs. The dvalues of the Obs are used as x- and yerror for the fit.
func (object): func has to be of the form
import autograd.numpy as anp def func(a, x): return a[0] + a[1] * x + a[2] * anp.sinh(x)
For multiple x values func can be of the form
def func(a, x): (x1, x2) = x return a[0] * x1 ** 2 + a[1] * x2
It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation will not work.
- silent (bool, optional): If true all output to the console is omitted (default False).
- initial_guess (list): can provide an initial guess for the input parameters. Relevant for non-linear fits with many parameters.
- expected_chisquare (bool): If true prints the expected chisquare which is corrected by effects caused by correlated input data. This can take a while as the full correlation matrix has to be calculated (default False).
Notes
Based on the orthogonal distance regression module of scipy
605def fit_lin(x, y, **kwargs): 606 """Performs a linear fit to y = n + m * x and returns two Obs n, m. 607 608 Parameters 609 ---------- 610 x : list 611 Can either be a list of floats in which case no xerror is assumed, or 612 a list of Obs, where the dvalues of the Obs are used as xerror for the fit. 613 y : list 614 List of Obs, the dvalues of the Obs are used as yerror for the fit. 615 """ 616 617 def f(a, x): 618 y = a[0] + a[1] * x 619 return y 620 621 if all(isinstance(n, Obs) for n in x): 622 out = total_least_squares(x, y, f, **kwargs) 623 return out.fit_parameters 624 elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray): 625 out = least_squares(x, y, f, **kwargs) 626 return out.fit_parameters 627 else: 628 raise Exception('Unsupported types for x')
Performs a linear fit to y = n + m * x and returns two Obs n, m.
Parameters
- x (list): Can either be a list of floats in which case no xerror is assumed, or a list of Obs, where the dvalues of the Obs are used as xerror for the fit.
- y (list): List of Obs, the dvalues of the Obs are used as yerror for the fit.
631def qqplot(x, o_y, func, p): 632 """Generates a quantile-quantile plot of the fit result which can be used to 633 check if the residuals of the fit are gaussian distributed. 634 """ 635 636 residuals = [] 637 for i_x, i_y in zip(x, o_y): 638 residuals.append((i_y - func(p, i_x)) / i_y.dvalue) 639 residuals = sorted(residuals) 640 my_y = [o.value for o in residuals] 641 probplot = scipy.stats.probplot(my_y) 642 my_x = probplot[0][0] 643 plt.figure(figsize=(8, 8 / 1.618)) 644 plt.errorbar(my_x, my_y, fmt='o') 645 fit_start = my_x[0] 646 fit_stop = my_x[-1] 647 samples = np.arange(fit_start, fit_stop, 0.01) 648 plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution') 649 plt.plot(samples, probplot[1][0] * samples + probplot[1][1], zorder=10, label='Least squares fit, r=' + str(np.around(probplot[1][2], 3)), marker='', ls='-') 650 651 plt.xlabel('Theoretical quantiles') 652 plt.ylabel('Ordered Values') 653 plt.legend() 654 plt.draw()
Generates a quantile-quantile plot of the fit result which can be used to check if the residuals of the fit are gaussian distributed.
657def residual_plot(x, y, func, fit_res): 658 """ Generates a plot which compares the fit to the data and displays the corresponding residuals""" 659 sorted_x = sorted(x) 660 xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0]) 661 xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2]) 662 x_samples = np.arange(xstart, xstop + 0.01, 0.01) 663 664 plt.figure(figsize=(8, 8 / 1.618)) 665 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0) 666 ax0 = plt.subplot(gs[0]) 667 ax0.errorbar(x, [o.value for o in y], yerr=[o.dvalue for o in y], ls='none', fmt='o', capsize=3, markersize=5, label='Data') 668 ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0) 669 ax0.set_xticklabels([]) 670 ax0.set_xlim([xstart, xstop]) 671 ax0.set_xticklabels([]) 672 ax0.legend() 673 674 residuals = (np.asarray([o.value for o in y]) - func([o.value for o in fit_res], x)) / np.asarray([o.dvalue for o in y]) 675 ax1 = plt.subplot(gs[1]) 676 ax1.plot(x, residuals, 'ko', ls='none', markersize=5) 677 ax1.tick_params(direction='out') 678 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True) 679 ax1.axhline(y=0.0, ls='--', color='k', marker=" ") 680 ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k') 681 ax1.set_xlim([xstart, xstop]) 682 ax1.set_ylabel('Residuals') 683 plt.subplots_adjust(wspace=None, hspace=None) 684 plt.draw()
Generates a plot which compares the fit to the data and displays the corresponding residuals
687def error_band(x, func, beta): 688 """Returns the error band for an array of sample values x, for given fit function func with optimized parameters beta.""" 689 cov = covariance(beta) 690 if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps): 691 warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning) 692 693 deriv = [] 694 for i, item in enumerate(x): 695 deriv.append(np.array(egrad(func)([o.value for o in beta], item))) 696 697 err = [] 698 for i, item in enumerate(x): 699 err.append(np.sqrt(deriv[i] @ cov @ deriv[i])) 700 err = np.array(err) 701 702 return err
Returns the error band for an array of sample values x, for given fit function func with optimized parameters beta.
705def ks_test(objects=None): 706 """Performs a Kolmogorov–Smirnov test for the p-values of all fit object. 707 708 Parameters 709 ---------- 710 objects : list 711 List of fit results to include in the analysis (optional). 712 """ 713 714 if objects is None: 715 obs_list = [] 716 for obj in gc.get_objects(): 717 if isinstance(obj, Fit_result): 718 obs_list.append(obj) 719 else: 720 obs_list = objects 721 722 p_values = [o.p_value for o in obs_list] 723 724 bins = len(p_values) 725 x = np.arange(0, 1.001, 0.001) 726 plt.plot(x, x, 'k', zorder=1) 727 plt.xlim(0, 1) 728 plt.ylim(0, 1) 729 plt.xlabel('p-value') 730 plt.ylabel('Cumulative probability') 731 plt.title(str(bins) + ' p-values') 732 733 n = np.arange(1, bins + 1) / np.float64(bins) 734 Xs = np.sort(p_values) 735 plt.step(Xs, n) 736 diffs = n - Xs 737 loc_max_diff = np.argmax(np.abs(diffs)) 738 loc = Xs[loc_max_diff] 739 plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0)) 740 plt.draw() 741 742 print(scipy.stats.kstest(p_values, 'uniform'))
Performs a Kolmogorov–Smirnov test for the p-values of all fit object.
Parameters
- objects (list): List of fit results to include in the analysis (optional).