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 as auto_jacobian 13from autograd import hessian as auto_hessian 14from autograd import elementwise_grad as egrad 15from numdifftools import Jacobian as num_jacobian 16from numdifftools import Hessian as num_hessian 17from .obs import Obs, derived_observable, covariance, cov_Obs, invert_corr_cov_cholesky 18 19 20class Fit_result(Sequence): 21 """Represents fit results. 22 23 Attributes 24 ---------- 25 fit_parameters : list 26 results for the individual fit parameters, 27 also accessible via indices. 28 chisquare_by_dof : float 29 reduced chisquare. 30 p_value : float 31 p-value of the fit 32 t2_p_value : float 33 Hotelling t-squared p-value for correlated fits. 34 """ 35 36 def __init__(self): 37 self.fit_parameters = None 38 39 def __getitem__(self, idx): 40 return self.fit_parameters[idx] 41 42 def __len__(self): 43 return len(self.fit_parameters) 44 45 def gamma_method(self, **kwargs): 46 """Apply the gamma method to all fit parameters""" 47 [o.gamma_method(**kwargs) for o in self.fit_parameters] 48 49 gm = gamma_method 50 51 def __str__(self): 52 my_str = 'Goodness of fit:\n' 53 if hasattr(self, 'chisquare_by_dof'): 54 my_str += '\u03C7\u00b2/d.o.f. = ' + f'{self.chisquare_by_dof:2.6f}' + '\n' 55 elif hasattr(self, 'residual_variance'): 56 my_str += 'residual variance = ' + f'{self.residual_variance:2.6f}' + '\n' 57 if hasattr(self, 'chisquare_by_expected_chisquare'): 58 my_str += '\u03C7\u00b2/\u03C7\u00b2exp = ' + f'{self.chisquare_by_expected_chisquare:2.6f}' + '\n' 59 if hasattr(self, 'p_value'): 60 my_str += 'p-value = ' + f'{self.p_value:2.4f}' + '\n' 61 if hasattr(self, 't2_p_value'): 62 my_str += 't\u00B2p-value = ' + f'{self.t2_p_value:2.4f}' + '\n' 63 my_str += 'Fit parameters:\n' 64 for i_par, par in enumerate(self.fit_parameters): 65 my_str += str(i_par) + '\t' + ' ' * int(par >= 0) + str(par).rjust(int(par < 0.0)) + '\n' 66 return my_str 67 68 def __repr__(self): 69 m = max(map(len, list(self.__dict__.keys()))) + 1 70 return '\n'.join([key.rjust(m) + ': ' + repr(value) for key, value in sorted(self.__dict__.items())]) 71 72 73def least_squares(x, y, func, priors=None, silent=False, **kwargs): 74 r'''Performs a non-linear fit to y = func(x). 75 ``` 76 77 Parameters 78 ---------- 79 For an uncombined fit: 80 81 x : list 82 list of floats. 83 y : list 84 list of Obs. 85 func : object 86 fit function, has to be of the form 87 88 ```python 89 import autograd.numpy as anp 90 91 def func(a, x): 92 return a[0] + a[1] * x + a[2] * anp.sinh(x) 93 ``` 94 95 For multiple x values func can be of the form 96 97 ```python 98 def func(a, x): 99 (x1, x2) = x 100 return a[0] * x1 ** 2 + a[1] * x2 101 ``` 102 It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation 103 will not work. 104 105 OR For a combined fit: 106 107 x : dict 108 dict of lists. 109 y : dict 110 dict of lists of Obs. 111 funcs : dict 112 dict of objects 113 fit functions have to be of the form (here a[0] is the common fit parameter) 114 ```python 115 import autograd.numpy as anp 116 funcs = {"a": func_a, 117 "b": func_b} 118 119 def func_a(a, x): 120 return a[1] * anp.exp(-a[0] * x) 121 122 def func_b(a, x): 123 return a[2] * anp.exp(-a[0] * x) 124 125 It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation 126 will not work. 127 128 priors : dict or list, optional 129 priors can either be a dictionary with integer keys and the corresponding priors as values or 130 a list with an entry for every parameter in the fit. The entries can either be 131 Obs (e.g. results from a previous fit) or strings containing a value and an error formatted like 132 0.548(23), 500(40) or 0.5(0.4) 133 silent : bool, optional 134 If true all output to the console is omitted (default False). 135 initial_guess : list 136 can provide an initial guess for the input parameters. Relevant for 137 non-linear fits with many parameters. In case of correlated fits the guess is used to perform 138 an uncorrelated fit which then serves as guess for the correlated fit. 139 method : str, optional 140 can be used to choose an alternative method for the minimization of chisquare. 141 The possible methods are the ones which can be used for scipy.optimize.minimize and 142 migrad of iminuit. If no method is specified, Levenberg-Marquard is used. 143 Reliable alternatives are migrad, Powell and Nelder-Mead. 144 tol: float, optional 145 can be used (only for combined fits and methods other than Levenberg-Marquard) to set the tolerance for convergence 146 to a different value to either speed up convergence at the cost of a larger error on the fitted parameters (and possibly 147 invalid estimates for parameter uncertainties) or smaller values to get more accurate parameter values 148 The stopping criterion depends on the method, e.g. migrad: edm_max = 0.002 * tol * errordef (EDM criterion: edm < edm_max) 149 correlated_fit : bool 150 If True, use the full inverse covariance matrix in the definition of the chisquare cost function. 151 For details about how the covariance matrix is estimated see `pyerrors.obs.covariance`. 152 In practice the correlation matrix is Cholesky decomposed and inverted (instead of the covariance matrix). 153 This procedure should be numerically more stable as the correlation matrix is typically better conditioned (Jacobi preconditioning). 154 inv_chol_cov_matrix [array,list], optional 155 array: shape = (no of y values) X (no of y values) 156 list: for an uncombined fit: [""] 157 for a combined fit: list of keys belonging to the corr_matrix saved in the array, must be the same as the keys of the y dict in alphabetical order 158 If correlated_fit=True is set as well, can provide an inverse covariance matrix (y errors, dy_f included!) of your own choosing for a correlated fit. 159 The matrix must be a lower triangular matrix constructed from a Cholesky decomposition: The function invert_corr_cov_cholesky(corr, inverrdiag) can be 160 used to construct it from a correlation matrix (corr) and the errors dy_f of the data points (inverrdiag = np.diag(1 / np.asarray(dy_f))). For the correct 161 ordering the correlation matrix (corr) can be sorted via the function sort_corr(corr, kl, yd) where kl is the list of keys and yd the y dict. 162 expected_chisquare : bool 163 If True estimates the expected chisquare which is 164 corrected by effects caused by correlated input data (default False). 165 resplot : bool 166 If True, a plot which displays fit, data and residuals is generated (default False). 167 qqplot : bool 168 If True, a quantile-quantile plot of the fit result is generated (default False). 169 num_grad : bool 170 Use numerical differentation instead of automatic differentiation to perform the error propagation (default False). 171 172 Returns 173 ------- 174 output : Fit_result 175 Parameters and information on the fitted result. 176 Examples 177 ------ 178 >>> # Example of a correlated (correlated_fit = True, inv_chol_cov_matrix handed over) combined fit, based on a randomly generated data set 179 >>> import numpy as np 180 >>> from scipy.stats import norm 181 >>> from scipy.linalg import cholesky 182 >>> import pyerrors as pe 183 >>> # generating the random data set 184 >>> num_samples = 400 185 >>> N = 3 186 >>> x = np.arange(N) 187 >>> x1 = norm.rvs(size=(N, num_samples)) # generate random numbers 188 >>> x2 = norm.rvs(size=(N, num_samples)) # generate random numbers 189 >>> r = r1 = r2 = np.zeros((N, N)) 190 >>> y = {} 191 >>> for i in range(N): 192 >>> for j in range(N): 193 >>> r[i, j] = np.exp(-0.8 * np.fabs(i - j)) # element in correlation matrix 194 >>> errl = np.sqrt([3.4, 2.5, 3.6]) # set y errors 195 >>> for i in range(N): 196 >>> for j in range(N): 197 >>> r[i, j] *= errl[i] * errl[j] # element in covariance matrix 198 >>> c = cholesky(r, lower=True) 199 >>> y = {'a': np.dot(c, x1), 'b': np.dot(c, x2)} # generate y data with the covariance matrix defined 200 >>> # random data set has been generated, now the dictionaries and the inverse covariance matrix to be handed over are built 201 >>> x_dict = {} 202 >>> y_dict = {} 203 >>> chol_inv_dict = {} 204 >>> data = [] 205 >>> for key in y.keys(): 206 >>> x_dict[key] = x 207 >>> for i in range(N): 208 >>> data.append(pe.Obs([[i + 1 + o for o in y[key][i]]], ['ens'])) # generate y Obs from the y data 209 >>> [o.gamma_method() for o in data] 210 >>> corr = pe.covariance(data, correlation=True) 211 >>> inverrdiag = np.diag(1 / np.asarray([o.dvalue for o in data])) 212 >>> chol_inv = pe.obs.invert_corr_cov_cholesky(corr, inverrdiag) # gives form of the inverse covariance matrix needed for the combined correlated fit below 213 >>> y_dict = {'a': data[:3], 'b': data[3:]} 214 >>> # common fit parameter p[0] in combined fit 215 >>> def fit1(p, x): 216 >>> return p[0] + p[1] * x 217 >>> def fit2(p, x): 218 >>> return p[0] + p[2] * x 219 >>> fitf_dict = {'a': fit1, 'b':fit2} 220 >>> fitp_inv_cov_combined_fit = pe.least_squares(x_dict,y_dict, fitf_dict, correlated_fit = True, inv_chol_cov_matrix = [chol_inv,['a','b']]) 221 Fit with 3 parameters 222 Method: Levenberg-Marquardt 223 `ftol` termination condition is satisfied. 224 chisquare/d.o.f.: 0.5388013574561786 # random 225 fit parameters [1.11897846 0.96361162 0.92325319] # random 226 227 ''' 228 output = Fit_result() 229 230 if (isinstance(x, dict) and isinstance(y, dict) and isinstance(func, dict)): 231 xd = {key: anp.asarray(x[key]) for key in x} 232 yd = y 233 funcd = func 234 output.fit_function = func 235 elif (isinstance(x, dict) or isinstance(y, dict) or isinstance(func, dict)): 236 raise TypeError("All arguments have to be dictionaries in order to perform a combined fit.") 237 else: 238 x = np.asarray(x) 239 xd = {"": x} 240 yd = {"": y} 241 funcd = {"": func} 242 output.fit_function = func 243 244 if kwargs.get('num_grad') is True: 245 jacobian = num_jacobian 246 hessian = num_hessian 247 else: 248 jacobian = auto_jacobian 249 hessian = auto_hessian 250 251 key_ls = sorted(list(xd.keys())) 252 253 if sorted(list(yd.keys())) != key_ls: 254 raise ValueError('x and y dictionaries do not contain the same keys.') 255 256 if sorted(list(funcd.keys())) != key_ls: 257 raise ValueError('x and func dictionaries do not contain the same keys.') 258 259 x_all = np.concatenate([np.array(xd[key]).transpose() for key in key_ls]).transpose() 260 y_all = np.concatenate([np.array(yd[key]) for key in key_ls]) 261 262 y_f = [o.value for o in y_all] 263 dy_f = [o.dvalue for o in y_all] 264 265 if len(x_all.shape) > 2: 266 raise ValueError("Unknown format for x values") 267 268 if np.any(np.asarray(dy_f) <= 0.0): 269 raise Exception("No y errors available, run the gamma method first.") 270 271 # number of fit parameters 272 n_parms_ls = [] 273 for key in key_ls: 274 if not callable(funcd[key]): 275 raise TypeError('func (key=' + key + ') is not a function.') 276 if np.asarray(xd[key]).shape[-1] != len(yd[key]): 277 raise ValueError('x and y input (key=' + key + ') do not have the same length') 278 for n_loc in range(100): 279 try: 280 funcd[key](np.arange(n_loc), x_all.T[0]) 281 except TypeError: 282 continue 283 except IndexError: 284 continue 285 else: 286 break 287 else: 288 raise RuntimeError("Fit function (key=" + key + ") is not valid.") 289 n_parms_ls.append(n_loc) 290 291 n_parms = max(n_parms_ls) 292 293 if len(key_ls) > 1: 294 for key in key_ls: 295 if np.asarray(yd[key]).shape != funcd[key](np.arange(n_parms), xd[key]).shape: 296 raise ValueError(f"Fit function {key} returns the wrong shape ({funcd[key](np.arange(n_parms), xd[key]).shape} instead of {xd[key].shape})\nIf the fit function is just a constant you could try adding x*0 to get the correct shape.") 297 298 if not silent: 299 print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1)) 300 301 if priors is not None: 302 if isinstance(priors, (list, np.ndarray)): 303 if n_parms != len(priors): 304 raise ValueError("'priors' does not have the correct length.") 305 306 loc_priors = [] 307 for i_n, i_prior in enumerate(priors): 308 loc_priors.append(_construct_prior_obs(i_prior, i_n)) 309 310 prior_mask = np.arange(len(priors)) 311 output.priors = loc_priors 312 313 elif isinstance(priors, dict): 314 loc_priors = [] 315 prior_mask = [] 316 output.priors = {} 317 for pos, prior in priors.items(): 318 if isinstance(pos, int): 319 prior_mask.append(pos) 320 else: 321 raise TypeError("Prior position needs to be an integer.") 322 loc_priors.append(_construct_prior_obs(prior, pos)) 323 324 output.priors[pos] = loc_priors[-1] 325 if max(prior_mask) >= n_parms: 326 raise ValueError("Prior position out of range.") 327 else: 328 raise TypeError("Unkown type for `priors`.") 329 330 p_f = [o.value for o in loc_priors] 331 dp_f = [o.dvalue for o in loc_priors] 332 if np.any(np.asarray(dp_f) <= 0.0): 333 raise Exception("No prior errors available, run the gamma method first.") 334 else: 335 p_f = dp_f = np.array([]) 336 prior_mask = [] 337 loc_priors = [] 338 339 if 'initial_guess' in kwargs: 340 x0 = kwargs.get('initial_guess') 341 if len(x0) != n_parms: 342 raise ValueError('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms)) 343 else: 344 x0 = [0.1] * n_parms 345 346 if priors is None: 347 def general_chisqfunc_uncorr(p, ivars, pr): 348 model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls]) 349 return (ivars - model) / dy_f 350 else: 351 def general_chisqfunc_uncorr(p, ivars, pr): 352 model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls]) 353 return anp.concatenate(((ivars - model) / dy_f, (p[prior_mask] - pr) / dp_f)) 354 355 def chisqfunc_uncorr(p): 356 return anp.sum(general_chisqfunc_uncorr(p, y_f, p_f) ** 2) 357 358 if kwargs.get('correlated_fit') is True: 359 if 'inv_chol_cov_matrix' in kwargs: 360 chol_inv = kwargs.get('inv_chol_cov_matrix') 361 if (chol_inv[0].shape[0] != len(dy_f)): 362 raise TypeError('The number of columns of the inverse covariance matrix handed over needs to be equal to the number of y errors.') 363 if (chol_inv[0].shape[0] != chol_inv[0].shape[1]): 364 raise TypeError('The inverse covariance matrix handed over needs to have the same number of rows as columns.') 365 if (chol_inv[1] != key_ls): 366 raise ValueError('The keys of inverse covariance matrix are not the same or do not appear in the same order as the x and y values.') 367 chol_inv = chol_inv[0] 368 else: 369 corr = covariance(y_all, correlation=True, **kwargs) 370 inverrdiag = np.diag(1 / np.asarray(dy_f)) 371 chol_inv = invert_corr_cov_cholesky(corr, inverrdiag) 372 373 def general_chisqfunc(p, ivars, pr): 374 model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls]) 375 return anp.concatenate((anp.dot(chol_inv, (ivars - model)), (p[prior_mask] - pr) / dp_f)) 376 377 def chisqfunc(p): 378 return anp.sum(general_chisqfunc(p, y_f, p_f) ** 2) 379 else: 380 general_chisqfunc = general_chisqfunc_uncorr 381 chisqfunc = chisqfunc_uncorr 382 383 output.method = kwargs.get('method', 'Levenberg-Marquardt') 384 if not silent: 385 print('Method:', output.method) 386 387 if output.method != 'Levenberg-Marquardt': 388 if output.method == 'migrad': 389 tolerance = 1e-4 # default value of 1e-1 set by iminuit can be problematic 390 if 'tol' in kwargs: 391 tolerance = kwargs.get('tol') 392 fit_result = iminuit.minimize(chisqfunc_uncorr, x0, tol=tolerance) # Stopping criterion 0.002 * tol * errordef 393 if kwargs.get('correlated_fit') is True: 394 fit_result = iminuit.minimize(chisqfunc, fit_result.x, tol=tolerance) 395 output.iterations = fit_result.nfev 396 else: 397 tolerance = 1e-12 398 if 'tol' in kwargs: 399 tolerance = kwargs.get('tol') 400 fit_result = scipy.optimize.minimize(chisqfunc_uncorr, x0, method=kwargs.get('method'), tol=tolerance) 401 if kwargs.get('correlated_fit') is True: 402 fit_result = scipy.optimize.minimize(chisqfunc, fit_result.x, method=kwargs.get('method'), tol=tolerance) 403 output.iterations = fit_result.nit 404 405 chisquare = fit_result.fun 406 407 else: 408 if 'tol' in kwargs: 409 print('tol cannot be set for Levenberg-Marquardt') 410 411 def chisqfunc_residuals_uncorr(p): 412 return general_chisqfunc_uncorr(p, y_f, p_f) 413 414 fit_result = scipy.optimize.least_squares(chisqfunc_residuals_uncorr, x0, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15) 415 if kwargs.get('correlated_fit') is True: 416 def chisqfunc_residuals(p): 417 return general_chisqfunc(p, y_f, p_f) 418 419 fit_result = scipy.optimize.least_squares(chisqfunc_residuals, fit_result.x, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15) 420 421 chisquare = np.sum(fit_result.fun ** 2) 422 assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14) 423 424 output.iterations = fit_result.nfev 425 426 if not fit_result.success: 427 raise Exception('The minimization procedure did not converge.') 428 429 output.chisquare = chisquare 430 output.dof = y_all.shape[-1] - n_parms + len(loc_priors) 431 output.p_value = 1 - scipy.stats.chi2.cdf(output.chisquare, output.dof) 432 if output.dof > 0: 433 output.chisquare_by_dof = output.chisquare / output.dof 434 else: 435 output.chisquare_by_dof = float('nan') 436 437 output.message = fit_result.message 438 if not silent: 439 print(fit_result.message) 440 print('chisquare/d.o.f.:', output.chisquare_by_dof) 441 print('fit parameters', fit_result.x) 442 443 def prepare_hat_matrix(): 444 hat_vector = [] 445 for key in key_ls: 446 if (len(xd[key]) != 0): 447 hat_vector.append(jacobian(funcd[key])(fit_result.x, xd[key])) 448 hat_vector = [item for sublist in hat_vector for item in sublist] 449 return hat_vector 450 451 if kwargs.get('expected_chisquare') is True: 452 if kwargs.get('correlated_fit') is not True: 453 W = np.diag(1 / np.asarray(dy_f)) 454 cov = covariance(y_all) 455 hat_vector = prepare_hat_matrix() 456 A = W @ hat_vector 457 P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T 458 expected_chisquare = np.trace((np.identity(y_all.shape[-1]) - P_phi) @ W @ cov @ W) 459 output.chisquare_by_expected_chisquare = output.chisquare / expected_chisquare 460 if not silent: 461 print('chisquare/expected_chisquare:', output.chisquare_by_expected_chisquare) 462 463 fitp = fit_result.x 464 465 try: 466 hess = hessian(chisqfunc)(fitp) 467 except TypeError: 468 raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None 469 470 len_y = len(y_f) 471 472 def chisqfunc_compact(d): 473 return anp.sum(general_chisqfunc(d[:n_parms], d[n_parms: n_parms + len_y], d[n_parms + len_y:]) ** 2) 474 475 jac_jac_y = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f, p_f))) 476 477 # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv 478 try: 479 deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms, n_parms:]) 480 except np.linalg.LinAlgError: 481 raise Exception("Cannot invert hessian matrix.") 482 483 result = [] 484 for i in range(n_parms): 485 result.append(derived_observable(lambda x_all, **kwargs: (x_all[0] + np.finfo(np.float64).eps) / (y_all[0].value + np.finfo(np.float64).eps) * fitp[i], list(y_all) + loc_priors, man_grad=list(deriv_y[i]))) 486 487 output.fit_parameters = result 488 489 # Hotelling t-squared p-value for correlated fits. 490 if kwargs.get('correlated_fit') is True: 491 n_cov = np.min(np.vectorize(lambda x_all: x_all.N)(y_all)) 492 output.t2_p_value = 1 - scipy.stats.f.cdf((n_cov - output.dof) / (output.dof * (n_cov - 1)) * output.chisquare, 493 output.dof, n_cov - output.dof) 494 495 if kwargs.get('resplot') is True: 496 for key in key_ls: 497 residual_plot(xd[key], yd[key], funcd[key], result, title=key) 498 499 if kwargs.get('qqplot') is True: 500 for key in key_ls: 501 qqplot(xd[key], yd[key], funcd[key], result, title=key) 502 503 return output 504 505 506def total_least_squares(x, y, func, silent=False, **kwargs): 507 r'''Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters. 508 509 Parameters 510 ---------- 511 x : list 512 list of Obs, or a tuple of lists of Obs 513 y : list 514 list of Obs. The dvalues of the Obs are used as x- and yerror for the fit. 515 func : object 516 func has to be of the form 517 518 ```python 519 import autograd.numpy as anp 520 521 def func(a, x): 522 return a[0] + a[1] * x + a[2] * anp.sinh(x) 523 ``` 524 525 For multiple x values func can be of the form 526 527 ```python 528 def func(a, x): 529 (x1, x2) = x 530 return a[0] * x1 ** 2 + a[1] * x2 531 ``` 532 533 It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation 534 will not work. 535 silent : bool, optional 536 If true all output to the console is omitted (default False). 537 initial_guess : list 538 can provide an initial guess for the input parameters. Relevant for non-linear 539 fits with many parameters. 540 expected_chisquare : bool 541 If true prints the expected chisquare which is 542 corrected by effects caused by correlated input data. 543 This can take a while as the full correlation matrix 544 has to be calculated (default False). 545 num_grad : bool 546 Use numerical differentation instead of automatic differentiation to perform the error propagation (default False). 547 548 Notes 549 ----- 550 Based on the orthogonal distance regression module of scipy. 551 552 Returns 553 ------- 554 output : Fit_result 555 Parameters and information on the fitted result. 556 ''' 557 558 output = Fit_result() 559 560 output.fit_function = func 561 562 x = np.array(x) 563 564 x_shape = x.shape 565 566 if kwargs.get('num_grad') is True: 567 jacobian = num_jacobian 568 hessian = num_hessian 569 else: 570 jacobian = auto_jacobian 571 hessian = auto_hessian 572 573 if not callable(func): 574 raise TypeError('func has to be a function.') 575 576 for i in range(42): 577 try: 578 func(np.arange(i), x.T[0]) 579 except TypeError: 580 continue 581 except IndexError: 582 continue 583 else: 584 break 585 else: 586 raise RuntimeError("Fit function is not valid.") 587 588 n_parms = i 589 if not silent: 590 print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1)) 591 592 x_f = np.vectorize(lambda o: o.value)(x) 593 dx_f = np.vectorize(lambda o: o.dvalue)(x) 594 y_f = np.array([o.value for o in y]) 595 dy_f = np.array([o.dvalue for o in y]) 596 597 if np.any(np.asarray(dx_f) <= 0.0): 598 raise Exception('No x errors available, run the gamma method first.') 599 600 if np.any(np.asarray(dy_f) <= 0.0): 601 raise Exception('No y errors available, run the gamma method first.') 602 603 if 'initial_guess' in kwargs: 604 x0 = kwargs.get('initial_guess') 605 if len(x0) != n_parms: 606 raise Exception('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms)) 607 else: 608 x0 = [1] * n_parms 609 610 data = RealData(x_f, y_f, sx=dx_f, sy=dy_f) 611 model = Model(func) 612 odr = ODR(data, model, x0, partol=np.finfo(np.float64).eps) 613 odr.set_job(fit_type=0, deriv=1) 614 out = odr.run() 615 616 output.residual_variance = out.res_var 617 618 output.method = 'ODR' 619 620 output.message = out.stopreason 621 622 output.xplus = out.xplus 623 624 if not silent: 625 print('Method: ODR') 626 print(*out.stopreason) 627 print('Residual variance:', output.residual_variance) 628 629 if out.info > 3: 630 raise Exception('The minimization procedure did not converge.') 631 632 m = x_f.size 633 634 def odr_chisquare(p): 635 model = func(p[:n_parms], p[n_parms:].reshape(x_shape)) 636 chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((x_f - p[n_parms:].reshape(x_shape)) / dx_f) ** 2) 637 return chisq 638 639 if kwargs.get('expected_chisquare') is True: 640 W = np.diag(1 / np.asarray(np.concatenate((dy_f.ravel(), dx_f.ravel())))) 641 642 if kwargs.get('covariance') is not None: 643 cov = kwargs.get('covariance') 644 else: 645 cov = covariance(np.concatenate((y, x.ravel()))) 646 647 number_of_x_parameters = int(m / x_f.shape[-1]) 648 649 old_jac = jacobian(func)(out.beta, out.xplus) 650 fused_row1 = np.concatenate((old_jac, np.concatenate((number_of_x_parameters * [np.zeros(old_jac.shape)]), axis=0))) 651 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]))) 652 new_jac = np.concatenate((fused_row1, fused_row2), axis=1) 653 654 A = W @ new_jac 655 P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T 656 expected_chisquare = np.trace((np.identity(P_phi.shape[0]) - P_phi) @ W @ cov @ W) 657 if expected_chisquare <= 0.0: 658 warnings.warn("Negative expected_chisquare.", RuntimeWarning) 659 expected_chisquare = np.abs(expected_chisquare) 660 output.chisquare_by_expected_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel()))) / expected_chisquare 661 if not silent: 662 print('chisquare/expected_chisquare:', 663 output.chisquare_by_expected_chisquare) 664 665 fitp = out.beta 666 try: 667 hess = hessian(odr_chisquare)(np.concatenate((fitp, out.xplus.ravel()))) 668 except TypeError: 669 raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None 670 671 def odr_chisquare_compact_x(d): 672 model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape)) 673 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) 674 return chisq 675 676 jac_jac_x = hessian(odr_chisquare_compact_x)(np.concatenate((fitp, out.xplus.ravel(), x_f.ravel()))) 677 678 # Compute hess^{-1} @ jac_jac_x[:n_parms + m, n_parms + m:] using LAPACK dgesv 679 try: 680 deriv_x = -scipy.linalg.solve(hess, jac_jac_x[:n_parms + m, n_parms + m:]) 681 except np.linalg.LinAlgError: 682 raise Exception("Cannot invert hessian matrix.") 683 684 def odr_chisquare_compact_y(d): 685 model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape)) 686 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) 687 return chisq 688 689 jac_jac_y = hessian(odr_chisquare_compact_y)(np.concatenate((fitp, out.xplus.ravel(), y_f))) 690 691 # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv 692 try: 693 deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms + m, n_parms + m:]) 694 except np.linalg.LinAlgError: 695 raise Exception("Cannot invert hessian matrix.") 696 697 result = [] 698 for i in range(n_parms): 699 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]))) 700 701 output.fit_parameters = result 702 703 output.odr_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel()))) 704 output.dof = x.shape[-1] - n_parms 705 output.p_value = 1 - scipy.stats.chi2.cdf(output.odr_chisquare, output.dof) 706 707 return output 708 709 710def fit_lin(x, y, **kwargs): 711 """Performs a linear fit to y = n + m * x and returns two Obs n, m. 712 713 Parameters 714 ---------- 715 x : list 716 Can either be a list of floats in which case no xerror is assumed, or 717 a list of Obs, where the dvalues of the Obs are used as xerror for the fit. 718 y : list 719 List of Obs, the dvalues of the Obs are used as yerror for the fit. 720 721 Returns 722 ------- 723 fit_parameters : list[Obs] 724 LIist of fitted observables. 725 """ 726 727 def f(a, x): 728 y = a[0] + a[1] * x 729 return y 730 731 if all(isinstance(n, Obs) for n in x): 732 out = total_least_squares(x, y, f, **kwargs) 733 return out.fit_parameters 734 elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray): 735 out = least_squares(x, y, f, **kwargs) 736 return out.fit_parameters 737 else: 738 raise TypeError('Unsupported types for x') 739 740 741def qqplot(x, o_y, func, p, title=""): 742 """Generates a quantile-quantile plot of the fit result which can be used to 743 check if the residuals of the fit are gaussian distributed. 744 745 Returns 746 ------- 747 None 748 """ 749 750 residuals = [] 751 for i_x, i_y in zip(x, o_y): 752 residuals.append((i_y - func(p, i_x)) / i_y.dvalue) 753 residuals = sorted(residuals) 754 my_y = [o.value for o in residuals] 755 probplot = scipy.stats.probplot(my_y) 756 my_x = probplot[0][0] 757 plt.figure(figsize=(8, 8 / 1.618)) 758 plt.errorbar(my_x, my_y, fmt='o') 759 fit_start = my_x[0] 760 fit_stop = my_x[-1] 761 samples = np.arange(fit_start, fit_stop, 0.01) 762 plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution') 763 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='-') 764 765 plt.xlabel('Theoretical quantiles') 766 plt.ylabel('Ordered Values') 767 plt.legend(title=title) 768 plt.draw() 769 770 771def residual_plot(x, y, func, fit_res, title=""): 772 """Generates a plot which compares the fit to the data and displays the corresponding residuals 773 774 For uncorrelated data the residuals are expected to be distributed ~N(0,1). 775 776 Returns 777 ------- 778 None 779 """ 780 sorted_x = sorted(x) 781 xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0]) 782 xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2]) 783 x_samples = np.arange(xstart, xstop + 0.01, 0.01) 784 785 plt.figure(figsize=(8, 8 / 1.618)) 786 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0) 787 ax0 = plt.subplot(gs[0]) 788 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') 789 ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0) 790 ax0.set_xticklabels([]) 791 ax0.set_xlim([xstart, xstop]) 792 ax0.set_xticklabels([]) 793 ax0.legend(title=title) 794 795 residuals = (np.asarray([o.value for o in y]) - func([o.value for o in fit_res], np.asarray(x))) / np.asarray([o.dvalue for o in y]) 796 ax1 = plt.subplot(gs[1]) 797 ax1.plot(x, residuals, 'ko', ls='none', markersize=5) 798 ax1.tick_params(direction='out') 799 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True) 800 ax1.axhline(y=0.0, ls='--', color='k', marker=" ") 801 ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k') 802 ax1.set_xlim([xstart, xstop]) 803 ax1.set_ylabel('Residuals') 804 plt.subplots_adjust(wspace=None, hspace=None) 805 plt.draw() 806 807 808def error_band(x, func, beta): 809 """Calculate the error band for an array of sample values x, for given fit function func with optimized parameters beta. 810 811 Returns 812 ------- 813 err : np.array(Obs) 814 Error band for an array of sample values x 815 """ 816 cov = covariance(beta) 817 if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps): 818 warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning) 819 820 deriv = [] 821 for i, item in enumerate(x): 822 deriv.append(np.array(egrad(func)([o.value for o in beta], item))) 823 824 err = [] 825 for i, item in enumerate(x): 826 err.append(np.sqrt(deriv[i] @ cov @ deriv[i])) 827 err = np.array(err) 828 829 return err 830 831 832def ks_test(objects=None): 833 """Performs a Kolmogorov–Smirnov test for the p-values of all fit object. 834 835 Parameters 836 ---------- 837 objects : list 838 List of fit results to include in the analysis (optional). 839 840 Returns 841 ------- 842 None 843 """ 844 845 if objects is None: 846 obs_list = [] 847 for obj in gc.get_objects(): 848 if isinstance(obj, Fit_result): 849 obs_list.append(obj) 850 else: 851 obs_list = objects 852 853 p_values = [o.p_value for o in obs_list] 854 855 bins = len(p_values) 856 x = np.arange(0, 1.001, 0.001) 857 plt.plot(x, x, 'k', zorder=1) 858 plt.xlim(0, 1) 859 plt.ylim(0, 1) 860 plt.xlabel('p-value') 861 plt.ylabel('Cumulative probability') 862 plt.title(str(bins) + ' p-values') 863 864 n = np.arange(1, bins + 1) / np.float64(bins) 865 Xs = np.sort(p_values) 866 plt.step(Xs, n) 867 diffs = n - Xs 868 loc_max_diff = np.argmax(np.abs(diffs)) 869 loc = Xs[loc_max_diff] 870 plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0)) 871 plt.draw() 872 873 print(scipy.stats.kstest(p_values, 'uniform')) 874 875 876def _extract_val_and_dval(string): 877 split_string = string.split('(') 878 if '.' in split_string[0] and '.' not in split_string[1][:-1]: 879 factor = 10 ** -len(split_string[0].partition('.')[2]) 880 else: 881 factor = 1 882 return float(split_string[0]), float(split_string[1][:-1]) * factor 883 884 885def _construct_prior_obs(i_prior, i_n): 886 if isinstance(i_prior, Obs): 887 return i_prior 888 elif isinstance(i_prior, str): 889 loc_val, loc_dval = _extract_val_and_dval(i_prior) 890 return cov_Obs(loc_val, loc_dval ** 2, '#prior' + str(i_n) + f"_{np.random.randint(2147483647):010d}") 891 else: 892 raise TypeError("Prior entries need to be 'Obs' or 'str'.")
21class Fit_result(Sequence): 22 """Represents fit results. 23 24 Attributes 25 ---------- 26 fit_parameters : list 27 results for the individual fit parameters, 28 also accessible via indices. 29 chisquare_by_dof : float 30 reduced chisquare. 31 p_value : float 32 p-value of the fit 33 t2_p_value : float 34 Hotelling t-squared p-value for correlated fits. 35 """ 36 37 def __init__(self): 38 self.fit_parameters = None 39 40 def __getitem__(self, idx): 41 return self.fit_parameters[idx] 42 43 def __len__(self): 44 return len(self.fit_parameters) 45 46 def gamma_method(self, **kwargs): 47 """Apply the gamma method to all fit parameters""" 48 [o.gamma_method(**kwargs) for o in self.fit_parameters] 49 50 gm = gamma_method 51 52 def __str__(self): 53 my_str = 'Goodness of fit:\n' 54 if hasattr(self, 'chisquare_by_dof'): 55 my_str += '\u03C7\u00b2/d.o.f. = ' + f'{self.chisquare_by_dof:2.6f}' + '\n' 56 elif hasattr(self, 'residual_variance'): 57 my_str += 'residual variance = ' + f'{self.residual_variance:2.6f}' + '\n' 58 if hasattr(self, 'chisquare_by_expected_chisquare'): 59 my_str += '\u03C7\u00b2/\u03C7\u00b2exp = ' + f'{self.chisquare_by_expected_chisquare:2.6f}' + '\n' 60 if hasattr(self, 'p_value'): 61 my_str += 'p-value = ' + f'{self.p_value:2.4f}' + '\n' 62 if hasattr(self, 't2_p_value'): 63 my_str += 't\u00B2p-value = ' + f'{self.t2_p_value:2.4f}' + '\n' 64 my_str += 'Fit parameters:\n' 65 for i_par, par in enumerate(self.fit_parameters): 66 my_str += str(i_par) + '\t' + ' ' * int(par >= 0) + str(par).rjust(int(par < 0.0)) + '\n' 67 return my_str 68 69 def __repr__(self): 70 m = max(map(len, list(self.__dict__.keys()))) + 1 71 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.
- chisquare_by_dof (float): reduced chisquare.
- p_value (float): p-value of the fit
- t2_p_value (float): Hotelling t-squared p-value for correlated fits.
74def least_squares(x, y, func, priors=None, silent=False, **kwargs): 75 r'''Performs a non-linear fit to y = func(x). 76 ``` 77 78 Parameters 79 ---------- 80 For an uncombined fit: 81 82 x : list 83 list of floats. 84 y : list 85 list of Obs. 86 func : object 87 fit function, has to be of the form 88 89 ```python 90 import autograd.numpy as anp 91 92 def func(a, x): 93 return a[0] + a[1] * x + a[2] * anp.sinh(x) 94 ``` 95 96 For multiple x values func can be of the form 97 98 ```python 99 def func(a, x): 100 (x1, x2) = x 101 return a[0] * x1 ** 2 + a[1] * x2 102 ``` 103 It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation 104 will not work. 105 106 OR For a combined fit: 107 108 x : dict 109 dict of lists. 110 y : dict 111 dict of lists of Obs. 112 funcs : dict 113 dict of objects 114 fit functions have to be of the form (here a[0] is the common fit parameter) 115 ```python 116 import autograd.numpy as anp 117 funcs = {"a": func_a, 118 "b": func_b} 119 120 def func_a(a, x): 121 return a[1] * anp.exp(-a[0] * x) 122 123 def func_b(a, x): 124 return a[2] * anp.exp(-a[0] * x) 125 126 It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation 127 will not work. 128 129 priors : dict or list, optional 130 priors can either be a dictionary with integer keys and the corresponding priors as values or 131 a list with an entry for every parameter in the fit. The entries can either be 132 Obs (e.g. results from a previous fit) or strings containing a value and an error formatted like 133 0.548(23), 500(40) or 0.5(0.4) 134 silent : bool, optional 135 If true all output to the console is omitted (default False). 136 initial_guess : list 137 can provide an initial guess for the input parameters. Relevant for 138 non-linear fits with many parameters. In case of correlated fits the guess is used to perform 139 an uncorrelated fit which then serves as guess for the correlated fit. 140 method : str, optional 141 can be used to choose an alternative method for the minimization of chisquare. 142 The possible methods are the ones which can be used for scipy.optimize.minimize and 143 migrad of iminuit. If no method is specified, Levenberg-Marquard is used. 144 Reliable alternatives are migrad, Powell and Nelder-Mead. 145 tol: float, optional 146 can be used (only for combined fits and methods other than Levenberg-Marquard) to set the tolerance for convergence 147 to a different value to either speed up convergence at the cost of a larger error on the fitted parameters (and possibly 148 invalid estimates for parameter uncertainties) or smaller values to get more accurate parameter values 149 The stopping criterion depends on the method, e.g. migrad: edm_max = 0.002 * tol * errordef (EDM criterion: edm < edm_max) 150 correlated_fit : bool 151 If True, use the full inverse covariance matrix in the definition of the chisquare cost function. 152 For details about how the covariance matrix is estimated see `pyerrors.obs.covariance`. 153 In practice the correlation matrix is Cholesky decomposed and inverted (instead of the covariance matrix). 154 This procedure should be numerically more stable as the correlation matrix is typically better conditioned (Jacobi preconditioning). 155 inv_chol_cov_matrix [array,list], optional 156 array: shape = (no of y values) X (no of y values) 157 list: for an uncombined fit: [""] 158 for a combined fit: list of keys belonging to the corr_matrix saved in the array, must be the same as the keys of the y dict in alphabetical order 159 If correlated_fit=True is set as well, can provide an inverse covariance matrix (y errors, dy_f included!) of your own choosing for a correlated fit. 160 The matrix must be a lower triangular matrix constructed from a Cholesky decomposition: The function invert_corr_cov_cholesky(corr, inverrdiag) can be 161 used to construct it from a correlation matrix (corr) and the errors dy_f of the data points (inverrdiag = np.diag(1 / np.asarray(dy_f))). For the correct 162 ordering the correlation matrix (corr) can be sorted via the function sort_corr(corr, kl, yd) where kl is the list of keys and yd the y dict. 163 expected_chisquare : bool 164 If True estimates the expected chisquare which is 165 corrected by effects caused by correlated input data (default False). 166 resplot : bool 167 If True, a plot which displays fit, data and residuals is generated (default False). 168 qqplot : bool 169 If True, a quantile-quantile plot of the fit result is generated (default False). 170 num_grad : bool 171 Use numerical differentation instead of automatic differentiation to perform the error propagation (default False). 172 173 Returns 174 ------- 175 output : Fit_result 176 Parameters and information on the fitted result. 177 Examples 178 ------ 179 >>> # Example of a correlated (correlated_fit = True, inv_chol_cov_matrix handed over) combined fit, based on a randomly generated data set 180 >>> import numpy as np 181 >>> from scipy.stats import norm 182 >>> from scipy.linalg import cholesky 183 >>> import pyerrors as pe 184 >>> # generating the random data set 185 >>> num_samples = 400 186 >>> N = 3 187 >>> x = np.arange(N) 188 >>> x1 = norm.rvs(size=(N, num_samples)) # generate random numbers 189 >>> x2 = norm.rvs(size=(N, num_samples)) # generate random numbers 190 >>> r = r1 = r2 = np.zeros((N, N)) 191 >>> y = {} 192 >>> for i in range(N): 193 >>> for j in range(N): 194 >>> r[i, j] = np.exp(-0.8 * np.fabs(i - j)) # element in correlation matrix 195 >>> errl = np.sqrt([3.4, 2.5, 3.6]) # set y errors 196 >>> for i in range(N): 197 >>> for j in range(N): 198 >>> r[i, j] *= errl[i] * errl[j] # element in covariance matrix 199 >>> c = cholesky(r, lower=True) 200 >>> y = {'a': np.dot(c, x1), 'b': np.dot(c, x2)} # generate y data with the covariance matrix defined 201 >>> # random data set has been generated, now the dictionaries and the inverse covariance matrix to be handed over are built 202 >>> x_dict = {} 203 >>> y_dict = {} 204 >>> chol_inv_dict = {} 205 >>> data = [] 206 >>> for key in y.keys(): 207 >>> x_dict[key] = x 208 >>> for i in range(N): 209 >>> data.append(pe.Obs([[i + 1 + o for o in y[key][i]]], ['ens'])) # generate y Obs from the y data 210 >>> [o.gamma_method() for o in data] 211 >>> corr = pe.covariance(data, correlation=True) 212 >>> inverrdiag = np.diag(1 / np.asarray([o.dvalue for o in data])) 213 >>> chol_inv = pe.obs.invert_corr_cov_cholesky(corr, inverrdiag) # gives form of the inverse covariance matrix needed for the combined correlated fit below 214 >>> y_dict = {'a': data[:3], 'b': data[3:]} 215 >>> # common fit parameter p[0] in combined fit 216 >>> def fit1(p, x): 217 >>> return p[0] + p[1] * x 218 >>> def fit2(p, x): 219 >>> return p[0] + p[2] * x 220 >>> fitf_dict = {'a': fit1, 'b':fit2} 221 >>> fitp_inv_cov_combined_fit = pe.least_squares(x_dict,y_dict, fitf_dict, correlated_fit = True, inv_chol_cov_matrix = [chol_inv,['a','b']]) 222 Fit with 3 parameters 223 Method: Levenberg-Marquardt 224 `ftol` termination condition is satisfied. 225 chisquare/d.o.f.: 0.5388013574561786 # random 226 fit parameters [1.11897846 0.96361162 0.92325319] # random 227 228 ''' 229 output = Fit_result() 230 231 if (isinstance(x, dict) and isinstance(y, dict) and isinstance(func, dict)): 232 xd = {key: anp.asarray(x[key]) for key in x} 233 yd = y 234 funcd = func 235 output.fit_function = func 236 elif (isinstance(x, dict) or isinstance(y, dict) or isinstance(func, dict)): 237 raise TypeError("All arguments have to be dictionaries in order to perform a combined fit.") 238 else: 239 x = np.asarray(x) 240 xd = {"": x} 241 yd = {"": y} 242 funcd = {"": func} 243 output.fit_function = func 244 245 if kwargs.get('num_grad') is True: 246 jacobian = num_jacobian 247 hessian = num_hessian 248 else: 249 jacobian = auto_jacobian 250 hessian = auto_hessian 251 252 key_ls = sorted(list(xd.keys())) 253 254 if sorted(list(yd.keys())) != key_ls: 255 raise ValueError('x and y dictionaries do not contain the same keys.') 256 257 if sorted(list(funcd.keys())) != key_ls: 258 raise ValueError('x and func dictionaries do not contain the same keys.') 259 260 x_all = np.concatenate([np.array(xd[key]).transpose() for key in key_ls]).transpose() 261 y_all = np.concatenate([np.array(yd[key]) for key in key_ls]) 262 263 y_f = [o.value for o in y_all] 264 dy_f = [o.dvalue for o in y_all] 265 266 if len(x_all.shape) > 2: 267 raise ValueError("Unknown format for x values") 268 269 if np.any(np.asarray(dy_f) <= 0.0): 270 raise Exception("No y errors available, run the gamma method first.") 271 272 # number of fit parameters 273 n_parms_ls = [] 274 for key in key_ls: 275 if not callable(funcd[key]): 276 raise TypeError('func (key=' + key + ') is not a function.') 277 if np.asarray(xd[key]).shape[-1] != len(yd[key]): 278 raise ValueError('x and y input (key=' + key + ') do not have the same length') 279 for n_loc in range(100): 280 try: 281 funcd[key](np.arange(n_loc), x_all.T[0]) 282 except TypeError: 283 continue 284 except IndexError: 285 continue 286 else: 287 break 288 else: 289 raise RuntimeError("Fit function (key=" + key + ") is not valid.") 290 n_parms_ls.append(n_loc) 291 292 n_parms = max(n_parms_ls) 293 294 if len(key_ls) > 1: 295 for key in key_ls: 296 if np.asarray(yd[key]).shape != funcd[key](np.arange(n_parms), xd[key]).shape: 297 raise ValueError(f"Fit function {key} returns the wrong shape ({funcd[key](np.arange(n_parms), xd[key]).shape} instead of {xd[key].shape})\nIf the fit function is just a constant you could try adding x*0 to get the correct shape.") 298 299 if not silent: 300 print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1)) 301 302 if priors is not None: 303 if isinstance(priors, (list, np.ndarray)): 304 if n_parms != len(priors): 305 raise ValueError("'priors' does not have the correct length.") 306 307 loc_priors = [] 308 for i_n, i_prior in enumerate(priors): 309 loc_priors.append(_construct_prior_obs(i_prior, i_n)) 310 311 prior_mask = np.arange(len(priors)) 312 output.priors = loc_priors 313 314 elif isinstance(priors, dict): 315 loc_priors = [] 316 prior_mask = [] 317 output.priors = {} 318 for pos, prior in priors.items(): 319 if isinstance(pos, int): 320 prior_mask.append(pos) 321 else: 322 raise TypeError("Prior position needs to be an integer.") 323 loc_priors.append(_construct_prior_obs(prior, pos)) 324 325 output.priors[pos] = loc_priors[-1] 326 if max(prior_mask) >= n_parms: 327 raise ValueError("Prior position out of range.") 328 else: 329 raise TypeError("Unkown type for `priors`.") 330 331 p_f = [o.value for o in loc_priors] 332 dp_f = [o.dvalue for o in loc_priors] 333 if np.any(np.asarray(dp_f) <= 0.0): 334 raise Exception("No prior errors available, run the gamma method first.") 335 else: 336 p_f = dp_f = np.array([]) 337 prior_mask = [] 338 loc_priors = [] 339 340 if 'initial_guess' in kwargs: 341 x0 = kwargs.get('initial_guess') 342 if len(x0) != n_parms: 343 raise ValueError('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms)) 344 else: 345 x0 = [0.1] * n_parms 346 347 if priors is None: 348 def general_chisqfunc_uncorr(p, ivars, pr): 349 model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls]) 350 return (ivars - model) / dy_f 351 else: 352 def general_chisqfunc_uncorr(p, ivars, pr): 353 model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls]) 354 return anp.concatenate(((ivars - model) / dy_f, (p[prior_mask] - pr) / dp_f)) 355 356 def chisqfunc_uncorr(p): 357 return anp.sum(general_chisqfunc_uncorr(p, y_f, p_f) ** 2) 358 359 if kwargs.get('correlated_fit') is True: 360 if 'inv_chol_cov_matrix' in kwargs: 361 chol_inv = kwargs.get('inv_chol_cov_matrix') 362 if (chol_inv[0].shape[0] != len(dy_f)): 363 raise TypeError('The number of columns of the inverse covariance matrix handed over needs to be equal to the number of y errors.') 364 if (chol_inv[0].shape[0] != chol_inv[0].shape[1]): 365 raise TypeError('The inverse covariance matrix handed over needs to have the same number of rows as columns.') 366 if (chol_inv[1] != key_ls): 367 raise ValueError('The keys of inverse covariance matrix are not the same or do not appear in the same order as the x and y values.') 368 chol_inv = chol_inv[0] 369 else: 370 corr = covariance(y_all, correlation=True, **kwargs) 371 inverrdiag = np.diag(1 / np.asarray(dy_f)) 372 chol_inv = invert_corr_cov_cholesky(corr, inverrdiag) 373 374 def general_chisqfunc(p, ivars, pr): 375 model = anp.concatenate([anp.array(funcd[key](p, xd[key])).reshape(-1) for key in key_ls]) 376 return anp.concatenate((anp.dot(chol_inv, (ivars - model)), (p[prior_mask] - pr) / dp_f)) 377 378 def chisqfunc(p): 379 return anp.sum(general_chisqfunc(p, y_f, p_f) ** 2) 380 else: 381 general_chisqfunc = general_chisqfunc_uncorr 382 chisqfunc = chisqfunc_uncorr 383 384 output.method = kwargs.get('method', 'Levenberg-Marquardt') 385 if not silent: 386 print('Method:', output.method) 387 388 if output.method != 'Levenberg-Marquardt': 389 if output.method == 'migrad': 390 tolerance = 1e-4 # default value of 1e-1 set by iminuit can be problematic 391 if 'tol' in kwargs: 392 tolerance = kwargs.get('tol') 393 fit_result = iminuit.minimize(chisqfunc_uncorr, x0, tol=tolerance) # Stopping criterion 0.002 * tol * errordef 394 if kwargs.get('correlated_fit') is True: 395 fit_result = iminuit.minimize(chisqfunc, fit_result.x, tol=tolerance) 396 output.iterations = fit_result.nfev 397 else: 398 tolerance = 1e-12 399 if 'tol' in kwargs: 400 tolerance = kwargs.get('tol') 401 fit_result = scipy.optimize.minimize(chisqfunc_uncorr, x0, method=kwargs.get('method'), tol=tolerance) 402 if kwargs.get('correlated_fit') is True: 403 fit_result = scipy.optimize.minimize(chisqfunc, fit_result.x, method=kwargs.get('method'), tol=tolerance) 404 output.iterations = fit_result.nit 405 406 chisquare = fit_result.fun 407 408 else: 409 if 'tol' in kwargs: 410 print('tol cannot be set for Levenberg-Marquardt') 411 412 def chisqfunc_residuals_uncorr(p): 413 return general_chisqfunc_uncorr(p, y_f, p_f) 414 415 fit_result = scipy.optimize.least_squares(chisqfunc_residuals_uncorr, x0, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15) 416 if kwargs.get('correlated_fit') is True: 417 def chisqfunc_residuals(p): 418 return general_chisqfunc(p, y_f, p_f) 419 420 fit_result = scipy.optimize.least_squares(chisqfunc_residuals, fit_result.x, method='lm', ftol=1e-15, gtol=1e-15, xtol=1e-15) 421 422 chisquare = np.sum(fit_result.fun ** 2) 423 assert np.isclose(chisquare, chisqfunc(fit_result.x), atol=1e-14) 424 425 output.iterations = fit_result.nfev 426 427 if not fit_result.success: 428 raise Exception('The minimization procedure did not converge.') 429 430 output.chisquare = chisquare 431 output.dof = y_all.shape[-1] - n_parms + len(loc_priors) 432 output.p_value = 1 - scipy.stats.chi2.cdf(output.chisquare, output.dof) 433 if output.dof > 0: 434 output.chisquare_by_dof = output.chisquare / output.dof 435 else: 436 output.chisquare_by_dof = float('nan') 437 438 output.message = fit_result.message 439 if not silent: 440 print(fit_result.message) 441 print('chisquare/d.o.f.:', output.chisquare_by_dof) 442 print('fit parameters', fit_result.x) 443 444 def prepare_hat_matrix(): 445 hat_vector = [] 446 for key in key_ls: 447 if (len(xd[key]) != 0): 448 hat_vector.append(jacobian(funcd[key])(fit_result.x, xd[key])) 449 hat_vector = [item for sublist in hat_vector for item in sublist] 450 return hat_vector 451 452 if kwargs.get('expected_chisquare') is True: 453 if kwargs.get('correlated_fit') is not True: 454 W = np.diag(1 / np.asarray(dy_f)) 455 cov = covariance(y_all) 456 hat_vector = prepare_hat_matrix() 457 A = W @ hat_vector 458 P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T 459 expected_chisquare = np.trace((np.identity(y_all.shape[-1]) - P_phi) @ W @ cov @ W) 460 output.chisquare_by_expected_chisquare = output.chisquare / expected_chisquare 461 if not silent: 462 print('chisquare/expected_chisquare:', output.chisquare_by_expected_chisquare) 463 464 fitp = fit_result.x 465 466 try: 467 hess = hessian(chisqfunc)(fitp) 468 except TypeError: 469 raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None 470 471 len_y = len(y_f) 472 473 def chisqfunc_compact(d): 474 return anp.sum(general_chisqfunc(d[:n_parms], d[n_parms: n_parms + len_y], d[n_parms + len_y:]) ** 2) 475 476 jac_jac_y = hessian(chisqfunc_compact)(np.concatenate((fitp, y_f, p_f))) 477 478 # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv 479 try: 480 deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms, n_parms:]) 481 except np.linalg.LinAlgError: 482 raise Exception("Cannot invert hessian matrix.") 483 484 result = [] 485 for i in range(n_parms): 486 result.append(derived_observable(lambda x_all, **kwargs: (x_all[0] + np.finfo(np.float64).eps) / (y_all[0].value + np.finfo(np.float64).eps) * fitp[i], list(y_all) + loc_priors, man_grad=list(deriv_y[i]))) 487 488 output.fit_parameters = result 489 490 # Hotelling t-squared p-value for correlated fits. 491 if kwargs.get('correlated_fit') is True: 492 n_cov = np.min(np.vectorize(lambda x_all: x_all.N)(y_all)) 493 output.t2_p_value = 1 - scipy.stats.f.cdf((n_cov - output.dof) / (output.dof * (n_cov - 1)) * output.chisquare, 494 output.dof, n_cov - output.dof) 495 496 if kwargs.get('resplot') is True: 497 for key in key_ls: 498 residual_plot(xd[key], yd[key], funcd[key], result, title=key) 499 500 if kwargs.get('qqplot') is True: 501 for key in key_ls: 502 qqplot(xd[key], yd[key], funcd[key], result, title=key) 503 504 return output
Performs a non-linear fit to y = func(x). ```
Parameters
- For an uncombined fit:
- 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.
- OR For a combined fit:
- x (dict): dict of lists.
- y (dict): dict of lists of Obs.
funcs (dict): dict of objects fit functions have to be of the form (here a[0] is the common fit parameter) ```python import autograd.numpy as anp funcs = {"a": func_a, "b": func_b}
def func_a(a, x): return a[1] * anp.exp(-a[0] * x)
def func_b(a, x): return a[2] * anp.exp(-a[0] * x)
It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation will not work.
- priors (dict or list, optional): priors can either be a dictionary with integer keys and the corresponding priors as values or 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.
- tol (float, optional): can be used (only for combined fits and methods other than Levenberg-Marquard) to set the tolerance for convergence to a different value to either speed up convergence at the cost of a larger error on the fitted parameters (and possibly invalid estimates for parameter uncertainties) or smaller values to get more accurate parameter values The stopping criterion depends on the method, e.g. migrad: edm_max = 0.002 * tol * errordef (EDM criterion: edm < edm_max)
- 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). - inv_chol_cov_matrix [array,list], optional: array: shape = (no of y values) X (no of y values) list: for an uncombined fit: [""] for a combined fit: list of keys belonging to the corr_matrix saved in the array, must be the same as the keys of the y dict in alphabetical order If correlated_fit=True is set as well, can provide an inverse covariance matrix (y errors, dy_f included!) of your own choosing for a correlated fit. The matrix must be a lower triangular matrix constructed from a Cholesky decomposition: The function invert_corr_cov_cholesky(corr, inverrdiag) can be used to construct it from a correlation matrix (corr) and the errors dy_f of the data points (inverrdiag = np.diag(1 / np.asarray(dy_f))). For the correct ordering the correlation matrix (corr) can be sorted via the function sort_corr(corr, kl, yd) where kl is the list of keys and yd the y dict.
- 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).
- num_grad (bool): Use numerical differentation instead of automatic differentiation to perform the error propagation (default False).
Returns
- output (Fit_result): Parameters and information on the fitted result.
Examples
>>> # Example of a correlated (correlated_fit = True, inv_chol_cov_matrix handed over) combined fit, based on a randomly generated data set
>>> import numpy as np
>>> from scipy.stats import norm
>>> from scipy.linalg import cholesky
>>> import pyerrors as pe
>>> # generating the random data set
>>> num_samples = 400
>>> N = 3
>>> x = np.arange(N)
>>> x1 = norm.rvs(size=(N, num_samples)) # generate random numbers
>>> x2 = norm.rvs(size=(N, num_samples)) # generate random numbers
>>> r = r1 = r2 = np.zeros((N, N))
>>> y = {}
>>> for i in range(N):
>>> for j in range(N):
>>> r[i, j] = np.exp(-0.8 * np.fabs(i - j)) # element in correlation matrix
>>> errl = np.sqrt([3.4, 2.5, 3.6]) # set y errors
>>> for i in range(N):
>>> for j in range(N):
>>> r[i, j] *= errl[i] * errl[j] # element in covariance matrix
>>> c = cholesky(r, lower=True)
>>> y = {'a': np.dot(c, x1), 'b': np.dot(c, x2)} # generate y data with the covariance matrix defined
>>> # random data set has been generated, now the dictionaries and the inverse covariance matrix to be handed over are built
>>> x_dict = {}
>>> y_dict = {}
>>> chol_inv_dict = {}
>>> data = []
>>> for key in y.keys():
>>> x_dict[key] = x
>>> for i in range(N):
>>> data.append(pe.Obs([[i + 1 + o for o in y[key][i]]], ['ens'])) # generate y Obs from the y data
>>> [o.gamma_method() for o in data]
>>> corr = pe.covariance(data, correlation=True)
>>> inverrdiag = np.diag(1 / np.asarray([o.dvalue for o in data]))
>>> chol_inv = pe.obs.invert_corr_cov_cholesky(corr, inverrdiag) # gives form of the inverse covariance matrix needed for the combined correlated fit below
>>> y_dict = {'a': data[:3], 'b': data[3:]}
>>> # common fit parameter p[0] in combined fit
>>> def fit1(p, x):
>>> return p[0] + p[1] * x
>>> def fit2(p, x):
>>> return p[0] + p[2] * x
>>> fitf_dict = {'a': fit1, 'b':fit2}
>>> fitp_inv_cov_combined_fit = pe.least_squares(x_dict,y_dict, fitf_dict, correlated_fit = True, inv_chol_cov_matrix = [chol_inv,['a','b']])
Fit with 3 parameters
Method: Levenberg-Marquardt
`ftol` termination condition is satisfied.
chisquare/d.o.f.: 0.5388013574561786 # random
fit parameters [1.11897846 0.96361162 0.92325319] # random
507def total_least_squares(x, y, func, silent=False, **kwargs): 508 r'''Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters. 509 510 Parameters 511 ---------- 512 x : list 513 list of Obs, or a tuple of lists of Obs 514 y : list 515 list of Obs. The dvalues of the Obs are used as x- and yerror for the fit. 516 func : object 517 func has to be of the form 518 519 ```python 520 import autograd.numpy as anp 521 522 def func(a, x): 523 return a[0] + a[1] * x + a[2] * anp.sinh(x) 524 ``` 525 526 For multiple x values func can be of the form 527 528 ```python 529 def func(a, x): 530 (x1, x2) = x 531 return a[0] * x1 ** 2 + a[1] * x2 532 ``` 533 534 It is important that all numpy functions refer to autograd.numpy, otherwise the differentiation 535 will not work. 536 silent : bool, optional 537 If true all output to the console is omitted (default False). 538 initial_guess : list 539 can provide an initial guess for the input parameters. Relevant for non-linear 540 fits with many parameters. 541 expected_chisquare : bool 542 If true prints the expected chisquare which is 543 corrected by effects caused by correlated input data. 544 This can take a while as the full correlation matrix 545 has to be calculated (default False). 546 num_grad : bool 547 Use numerical differentation instead of automatic differentiation to perform the error propagation (default False). 548 549 Notes 550 ----- 551 Based on the orthogonal distance regression module of scipy. 552 553 Returns 554 ------- 555 output : Fit_result 556 Parameters and information on the fitted result. 557 ''' 558 559 output = Fit_result() 560 561 output.fit_function = func 562 563 x = np.array(x) 564 565 x_shape = x.shape 566 567 if kwargs.get('num_grad') is True: 568 jacobian = num_jacobian 569 hessian = num_hessian 570 else: 571 jacobian = auto_jacobian 572 hessian = auto_hessian 573 574 if not callable(func): 575 raise TypeError('func has to be a function.') 576 577 for i in range(42): 578 try: 579 func(np.arange(i), x.T[0]) 580 except TypeError: 581 continue 582 except IndexError: 583 continue 584 else: 585 break 586 else: 587 raise RuntimeError("Fit function is not valid.") 588 589 n_parms = i 590 if not silent: 591 print('Fit with', n_parms, 'parameter' + 's' * (n_parms > 1)) 592 593 x_f = np.vectorize(lambda o: o.value)(x) 594 dx_f = np.vectorize(lambda o: o.dvalue)(x) 595 y_f = np.array([o.value for o in y]) 596 dy_f = np.array([o.dvalue for o in y]) 597 598 if np.any(np.asarray(dx_f) <= 0.0): 599 raise Exception('No x errors available, run the gamma method first.') 600 601 if np.any(np.asarray(dy_f) <= 0.0): 602 raise Exception('No y errors available, run the gamma method first.') 603 604 if 'initial_guess' in kwargs: 605 x0 = kwargs.get('initial_guess') 606 if len(x0) != n_parms: 607 raise Exception('Initial guess does not have the correct length: %d vs. %d' % (len(x0), n_parms)) 608 else: 609 x0 = [1] * n_parms 610 611 data = RealData(x_f, y_f, sx=dx_f, sy=dy_f) 612 model = Model(func) 613 odr = ODR(data, model, x0, partol=np.finfo(np.float64).eps) 614 odr.set_job(fit_type=0, deriv=1) 615 out = odr.run() 616 617 output.residual_variance = out.res_var 618 619 output.method = 'ODR' 620 621 output.message = out.stopreason 622 623 output.xplus = out.xplus 624 625 if not silent: 626 print('Method: ODR') 627 print(*out.stopreason) 628 print('Residual variance:', output.residual_variance) 629 630 if out.info > 3: 631 raise Exception('The minimization procedure did not converge.') 632 633 m = x_f.size 634 635 def odr_chisquare(p): 636 model = func(p[:n_parms], p[n_parms:].reshape(x_shape)) 637 chisq = anp.sum(((y_f - model) / dy_f) ** 2) + anp.sum(((x_f - p[n_parms:].reshape(x_shape)) / dx_f) ** 2) 638 return chisq 639 640 if kwargs.get('expected_chisquare') is True: 641 W = np.diag(1 / np.asarray(np.concatenate((dy_f.ravel(), dx_f.ravel())))) 642 643 if kwargs.get('covariance') is not None: 644 cov = kwargs.get('covariance') 645 else: 646 cov = covariance(np.concatenate((y, x.ravel()))) 647 648 number_of_x_parameters = int(m / x_f.shape[-1]) 649 650 old_jac = jacobian(func)(out.beta, out.xplus) 651 fused_row1 = np.concatenate((old_jac, np.concatenate((number_of_x_parameters * [np.zeros(old_jac.shape)]), axis=0))) 652 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]))) 653 new_jac = np.concatenate((fused_row1, fused_row2), axis=1) 654 655 A = W @ new_jac 656 P_phi = A @ np.linalg.pinv(A.T @ A) @ A.T 657 expected_chisquare = np.trace((np.identity(P_phi.shape[0]) - P_phi) @ W @ cov @ W) 658 if expected_chisquare <= 0.0: 659 warnings.warn("Negative expected_chisquare.", RuntimeWarning) 660 expected_chisquare = np.abs(expected_chisquare) 661 output.chisquare_by_expected_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel()))) / expected_chisquare 662 if not silent: 663 print('chisquare/expected_chisquare:', 664 output.chisquare_by_expected_chisquare) 665 666 fitp = out.beta 667 try: 668 hess = hessian(odr_chisquare)(np.concatenate((fitp, out.xplus.ravel()))) 669 except TypeError: 670 raise Exception("It is required to use autograd.numpy instead of numpy within fit functions, see the documentation for details.") from None 671 672 def odr_chisquare_compact_x(d): 673 model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape)) 674 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) 675 return chisq 676 677 jac_jac_x = hessian(odr_chisquare_compact_x)(np.concatenate((fitp, out.xplus.ravel(), x_f.ravel()))) 678 679 # Compute hess^{-1} @ jac_jac_x[:n_parms + m, n_parms + m:] using LAPACK dgesv 680 try: 681 deriv_x = -scipy.linalg.solve(hess, jac_jac_x[:n_parms + m, n_parms + m:]) 682 except np.linalg.LinAlgError: 683 raise Exception("Cannot invert hessian matrix.") 684 685 def odr_chisquare_compact_y(d): 686 model = func(d[:n_parms], d[n_parms:n_parms + m].reshape(x_shape)) 687 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) 688 return chisq 689 690 jac_jac_y = hessian(odr_chisquare_compact_y)(np.concatenate((fitp, out.xplus.ravel(), y_f))) 691 692 # Compute hess^{-1} @ jac_jac_y[:n_parms + m, n_parms + m:] using LAPACK dgesv 693 try: 694 deriv_y = -scipy.linalg.solve(hess, jac_jac_y[:n_parms + m, n_parms + m:]) 695 except np.linalg.LinAlgError: 696 raise Exception("Cannot invert hessian matrix.") 697 698 result = [] 699 for i in range(n_parms): 700 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]))) 701 702 output.fit_parameters = result 703 704 output.odr_chisquare = odr_chisquare(np.concatenate((out.beta, out.xplus.ravel()))) 705 output.dof = x.shape[-1] - n_parms 706 output.p_value = 1 - scipy.stats.chi2.cdf(output.odr_chisquare, output.dof) 707 708 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).
- num_grad (bool): Use numerical differentation instead of automatic differentiation to perform the error propagation (default False).
Notes
Based on the orthogonal distance regression module of scipy.
Returns
- output (Fit_result): Parameters and information on the fitted result.
711def fit_lin(x, y, **kwargs): 712 """Performs a linear fit to y = n + m * x and returns two Obs n, m. 713 714 Parameters 715 ---------- 716 x : list 717 Can either be a list of floats in which case no xerror is assumed, or 718 a list of Obs, where the dvalues of the Obs are used as xerror for the fit. 719 y : list 720 List of Obs, the dvalues of the Obs are used as yerror for the fit. 721 722 Returns 723 ------- 724 fit_parameters : list[Obs] 725 LIist of fitted observables. 726 """ 727 728 def f(a, x): 729 y = a[0] + a[1] * x 730 return y 731 732 if all(isinstance(n, Obs) for n in x): 733 out = total_least_squares(x, y, f, **kwargs) 734 return out.fit_parameters 735 elif all(isinstance(n, float) or isinstance(n, int) for n in x) or isinstance(x, np.ndarray): 736 out = least_squares(x, y, f, **kwargs) 737 return out.fit_parameters 738 else: 739 raise TypeError('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.
Returns
- fit_parameters (list[Obs]): LIist of fitted observables.
742def qqplot(x, o_y, func, p, title=""): 743 """Generates a quantile-quantile plot of the fit result which can be used to 744 check if the residuals of the fit are gaussian distributed. 745 746 Returns 747 ------- 748 None 749 """ 750 751 residuals = [] 752 for i_x, i_y in zip(x, o_y): 753 residuals.append((i_y - func(p, i_x)) / i_y.dvalue) 754 residuals = sorted(residuals) 755 my_y = [o.value for o in residuals] 756 probplot = scipy.stats.probplot(my_y) 757 my_x = probplot[0][0] 758 plt.figure(figsize=(8, 8 / 1.618)) 759 plt.errorbar(my_x, my_y, fmt='o') 760 fit_start = my_x[0] 761 fit_stop = my_x[-1] 762 samples = np.arange(fit_start, fit_stop, 0.01) 763 plt.plot(samples, samples, 'k--', zorder=11, label='Standard normal distribution') 764 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='-') 765 766 plt.xlabel('Theoretical quantiles') 767 plt.ylabel('Ordered Values') 768 plt.legend(title=title) 769 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.
Returns
- None
772def residual_plot(x, y, func, fit_res, title=""): 773 """Generates a plot which compares the fit to the data and displays the corresponding residuals 774 775 For uncorrelated data the residuals are expected to be distributed ~N(0,1). 776 777 Returns 778 ------- 779 None 780 """ 781 sorted_x = sorted(x) 782 xstart = sorted_x[0] - 0.5 * (sorted_x[1] - sorted_x[0]) 783 xstop = sorted_x[-1] + 0.5 * (sorted_x[-1] - sorted_x[-2]) 784 x_samples = np.arange(xstart, xstop + 0.01, 0.01) 785 786 plt.figure(figsize=(8, 8 / 1.618)) 787 gs = gridspec.GridSpec(2, 1, height_ratios=[3, 1], wspace=0.0, hspace=0.0) 788 ax0 = plt.subplot(gs[0]) 789 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') 790 ax0.plot(x_samples, func([o.value for o in fit_res], x_samples), label='Fit', zorder=10, ls='-', ms=0) 791 ax0.set_xticklabels([]) 792 ax0.set_xlim([xstart, xstop]) 793 ax0.set_xticklabels([]) 794 ax0.legend(title=title) 795 796 residuals = (np.asarray([o.value for o in y]) - func([o.value for o in fit_res], np.asarray(x))) / np.asarray([o.dvalue for o in y]) 797 ax1 = plt.subplot(gs[1]) 798 ax1.plot(x, residuals, 'ko', ls='none', markersize=5) 799 ax1.tick_params(direction='out') 800 ax1.tick_params(axis="x", bottom=True, top=True, labelbottom=True) 801 ax1.axhline(y=0.0, ls='--', color='k', marker=" ") 802 ax1.fill_between(x_samples, -1.0, 1.0, alpha=0.1, facecolor='k') 803 ax1.set_xlim([xstart, xstop]) 804 ax1.set_ylabel('Residuals') 805 plt.subplots_adjust(wspace=None, hspace=None) 806 plt.draw()
Generates a plot which compares the fit to the data and displays the corresponding residuals
For uncorrelated data the residuals are expected to be distributed ~N(0,1).
Returns
- None
809def error_band(x, func, beta): 810 """Calculate the error band for an array of sample values x, for given fit function func with optimized parameters beta. 811 812 Returns 813 ------- 814 err : np.array(Obs) 815 Error band for an array of sample values x 816 """ 817 cov = covariance(beta) 818 if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps): 819 warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning) 820 821 deriv = [] 822 for i, item in enumerate(x): 823 deriv.append(np.array(egrad(func)([o.value for o in beta], item))) 824 825 err = [] 826 for i, item in enumerate(x): 827 err.append(np.sqrt(deriv[i] @ cov @ deriv[i])) 828 err = np.array(err) 829 830 return err
Calculate the error band for an array of sample values x, for given fit function func with optimized parameters beta.
Returns
- err (np.array(Obs)): Error band for an array of sample values x
833def ks_test(objects=None): 834 """Performs a Kolmogorov–Smirnov test for the p-values of all fit object. 835 836 Parameters 837 ---------- 838 objects : list 839 List of fit results to include in the analysis (optional). 840 841 Returns 842 ------- 843 None 844 """ 845 846 if objects is None: 847 obs_list = [] 848 for obj in gc.get_objects(): 849 if isinstance(obj, Fit_result): 850 obs_list.append(obj) 851 else: 852 obs_list = objects 853 854 p_values = [o.p_value for o in obs_list] 855 856 bins = len(p_values) 857 x = np.arange(0, 1.001, 0.001) 858 plt.plot(x, x, 'k', zorder=1) 859 plt.xlim(0, 1) 860 plt.ylim(0, 1) 861 plt.xlabel('p-value') 862 plt.ylabel('Cumulative probability') 863 plt.title(str(bins) + ' p-values') 864 865 n = np.arange(1, bins + 1) / np.float64(bins) 866 Xs = np.sort(p_values) 867 plt.step(Xs, n) 868 diffs = n - Xs 869 loc_max_diff = np.argmax(np.abs(diffs)) 870 loc = Xs[loc_max_diff] 871 plt.annotate('', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0)) 872 plt.draw() 873 874 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).
Returns
- None