diff --git a/README.md b/README.md index 64efa51d..e67157ea 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,4 @@ -[![flake8 Lint](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml) [![CI](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml) [![](https://img.shields.io/badge/python-3.6+-blue.svg)](https://www.python.org/downloads/) +[![flake8](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/flake8.yml) [![CI](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml/badge.svg)](https://github.com/fjosw/pyerrors/actions/workflows/CI.yml) [![](https://img.shields.io/badge/python-3.6+-blue.svg)](https://www.python.org/downloads/) # pyerrors `pyerrors` is a python package for error computation and propagation of Markov chain Monte Carlo data. It is based on the **gamma method** [arXiv:hep-lat/0306017](https://arxiv.org/abs/hep-lat/0306017). Some of its features are: diff --git a/pyerrors/fits.py b/pyerrors/fits.py index d57445bd..fbaca972 100644 --- a/pyerrors/fits.py +++ b/pyerrors/fits.py @@ -1,6 +1,7 @@ #!/usr/bin/env python # coding: utf-8 +import gc import warnings import numpy as np import autograd.numpy as anp @@ -598,6 +599,52 @@ def error_band(x, func, beta): return err +def ks_test(obs=None): + """Performs a Kolmogorov–Smirnov test for the Q-values of all fit object. + + If no list is given all Obs in memory are used. + + Disclaimer: The determination of the individual Q-values as well as this function have not been tested yet. + """ + + raise Exception('Not yet implemented') + + if obs is None: + obs_list = [] + for obj in gc.get_objects(): + if isinstance(obj, Obs): + obs_list.append(obj) + else: + obs_list = obs + + # TODO: Rework to apply to Q-values of all fits in memory + Qs = [] + for obs_i in obs_list: + for ens in obs_i.e_names: + if obs_i.e_Q[ens] is not None: + Qs.append(obs_i.e_Q[ens]) + + bins = len(Qs) + x = np.arange(0, 1.001, 0.001) + plt.plot(x, x, 'k', zorder=1) + plt.xlim(0, 1) + plt.ylim(0, 1) + plt.xlabel('Q value') + plt.ylabel('Cumulative probability') + plt.title(str(bins) + ' Q values') + + n = np.arange(1, bins + 1) / np.float64(bins) + Xs = np.sort(Qs) + plt.step(Xs, n) + diffs = n - Xs + loc_max_diff = np.argmax(np.abs(diffs)) + loc = Xs[loc_max_diff] + plt.annotate(s='', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0)) + plt.show() + + print(scipy.stats.kstest(Qs, 'uniform')) + + def fit_general(x, y, func, silent=False, **kwargs): """Performs a non-linear fit to y = func(x) and returns a list of Obs corresponding to the fit parameters. diff --git a/pyerrors/misc.py b/pyerrors/misc.py index f122fac2..920518e5 100644 --- a/pyerrors/misc.py +++ b/pyerrors/misc.py @@ -1,10 +1,7 @@ #!/usr/bin/env python # coding: utf-8 -import gc import numpy as np -import scipy.stats -import matplotlib.pyplot as plt from .pyerrors import Obs @@ -38,47 +35,3 @@ def gen_correlated_data(means, cov, name, tau=0.5, samples=1000): data.append(np.sqrt(1 - a ** 2) * rand[i] + a * data[-1]) corr_data = np.array(data) - np.mean(data, axis=0) + means return [Obs([dat], [name]) for dat in corr_data.T] - - -def ks_test(obs=None): - """Performs a Kolmogorov–Smirnov test for the Q-values of a list of Obs. - - If no list is given all Obs in memory are used. - - Disclaimer: The determination of the individual Q-values as well as this function have not been tested yet. - """ - - if obs is None: - obs_list = [] - for obj in gc.get_objects(): - if isinstance(obj, Obs): - obs_list.append(obj) - else: - obs_list = obs - - # TODO: Rework to apply to Q-values of all fits in memory - Qs = [] - for obs_i in obs_list: - for ens in obs_i.e_names: - if obs_i.e_Q[ens] is not None: - Qs.append(obs_i.e_Q[ens]) - - bins = len(Qs) - x = np.arange(0, 1.001, 0.001) - plt.plot(x, x, 'k', zorder=1) - plt.xlim(0, 1) - plt.ylim(0, 1) - plt.xlabel('Q value') - plt.ylabel('Cumulative probability') - plt.title(str(bins) + ' Q values') - - n = np.arange(1, bins + 1) / np.float64(bins) - Xs = np.sort(Qs) - plt.step(Xs, n) - diffs = n - Xs - loc_max_diff = np.argmax(np.abs(diffs)) - loc = Xs[loc_max_diff] - plt.annotate(s='', xy=(loc, loc), xytext=(loc, loc + diffs[loc_max_diff]), arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0)) - plt.show() - - print(scipy.stats.kstest(Qs, 'uniform'))