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