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