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