From 0c60032924bd5d6b6aacade165d7074178cc2652 Mon Sep 17 00:00:00 2001 From: Patricio Cubillos Date: Wed, 16 Oct 2019 16:24:00 +0200 Subject: [PATCH 1/6] Fixed bug in stats.dwt_chisq() when including priors. Updated docstrings and added tests. --- mc3/stats/stats.py | 16 +++++++++++++--- tests/test_stats.py | 24 ++++++++++++++++++++++-- 2 files changed, 35 insertions(+), 5 deletions(-) diff --git a/mc3/stats/stats.py b/mc3/stats/stats.py index 35990dd..66b8cae 100644 --- a/mc3/stats/stats.py +++ b/mc3/stats/stats.py @@ -179,6 +179,7 @@ def chisq(model, data, uncert, Examples -------- >>> import mc3.stats as ms + >>> import numpy as np >>> # Compute chi-squared for a given model fitting a data set: >>> data = np.array([1.1, 1.2, 0.9, 1.0]) >>> model = np.array([1.0, 1.0, 1.0, 1.0]) @@ -243,20 +244,29 @@ def dwt_chisq(model, data, params, priors=None, priorlow=None, priorup=None): -------- >>> import mc3.stats as ms >>> import numpy as np - + >>> # Compute chi-squared for a given model fitting a data set: >>> data = np.array([2.0, 0.0, 3.0, -2.0, -1.0, 2.0, 2.0, 0.0]) >>> model = np.ones(8) >>> params = np.array([1.0, 0.1, 0.1]) - >>> chisq = ms.chisq(model, data, params) + >>> chisq = ms.dwt_chisq(model, data, params) >>> print(chisq) 1693.22308882 + >>> # Now, say this is a three-parameter model, with a Gaussian prior + >>> # on the last parameter: + >>> priors = np.array([1.0, 0.2, 0.3]) + >>> plow = np.array([0.0, 0.0, 0.1]) + >>> pup = np.array([0.0, 0.0, 0.1]) + >>> chisq = ms.dwt_chisq(model, data, params, priors, plow, pup) + >>> print(chisq) + 1697.2230888243134 """ if priors is None or priorlow is None or priorup is None: return dwt.chisq(params, model, data) iprior = (priorlow > 0) & (priorup > 0) dprior = (params - priors)[iprior] - return dwt.chisq(params, model, data, dprior, priorlow, priorup) + return dwt.chisq(params, model, data, dprior, + priorlow[iprior], priorup[iprior]) def log_prior(posterior, prior, priorlow, priorup, pstep): diff --git a/tests/test_stats.py b/tests/test_stats.py index cf17b93..94d340f 100644 --- a/tests/test_stats.py +++ b/tests/test_stats.py @@ -107,7 +107,7 @@ def test_residuals(): np.testing.assert_allclose(res, np.array([-1.0, -2.0, 1.0, 0.0, 0.5])) -def test_chisq_no_priors(): +def test_chisq(): data = np.array([1.1, 1.2, 0.9, 1.0]) model = np.array([1.0, 1.0, 1.0, 1.0]) uncert = np.array([0.1, 0.1, 0.1, 0.1]) @@ -115,7 +115,7 @@ def test_chisq_no_priors(): assert chisq == 6.0 -def test_chisq(): +def test_chisq_priors(): data = np.array([1.1, 1.2, 0.9, 1.0]) model = np.array([1.0, 1.0, 1.0, 1.0]) uncert = np.array([0.1, 0.1, 0.1, 0.1]) @@ -127,6 +127,26 @@ def test_chisq(): assert chisq == 6.25 +def test_dwt_chisq(): + data = np.array([2.0, 0.0, 3.0, -2.0, -1.0, 2.0, 2.0, 0.0]) + model = np.ones(8) + params = np.array([1.0, 0.1, 0.1]) + chisq = ms.dwt_chisq(model, data, params) + np.testing.assert_allclose(chisq, 1693.22308882) + + +def test_dwt_chisq_priors(): + data = np.array([2.0, 0.0, 3.0, -2.0, -1.0, 2.0, 2.0, 0.0]) + model = np.ones(8) + params = np.array([1.0, 0.1, 0.1]) + priors = np.array([1.0, 0.2, 0.3]) + plow = np.array([0.0, 0.0, 0.1]) + pup = np.array([0.0, 0.0, 0.1]) + chisq = ms.dwt_chisq(model, data, params, priors, plow, pup) + np.testing.assert_allclose(chisq, 1697.2230888243134) + + + def test_log_prior_uniform(): post = np.array([[3.0, 2.0], [3.1, 1.0], [3.6, 1.5]]) prior = np.array([3.5, 0.0]) From 2ed90af8a678147aa3b94a4a83af3b81e0de43a9 Mon Sep 17 00:00:00 2001 From: Patricio Cubillos Date: Wed, 16 Oct 2019 16:34:00 +0200 Subject: [PATCH 2/6] Raise error when stats.dwt_chisq() has an input params with less than two values, added test. --- mc3/stats/stats.py | 4 ++++ tests/test_stats.py | 7 +++++++ 2 files changed, 11 insertions(+) diff --git a/mc3/stats/stats.py b/mc3/stats/stats.py index 66b8cae..c872716 100644 --- a/mc3/stats/stats.py +++ b/mc3/stats/stats.py @@ -260,6 +260,10 @@ def dwt_chisq(model, data, params, priors=None, priorlow=None, priorup=None): >>> print(chisq) 1697.2230888243134 """ + if len(params) < 3: + with mu.Log() as log: + log.error('Wavelet chisq should have at least three parameters.') + if priors is None or priorlow is None or priorup is None: return dwt.chisq(params, model, data) diff --git a/tests/test_stats.py b/tests/test_stats.py index 94d340f..bb39323 100644 --- a/tests/test_stats.py +++ b/tests/test_stats.py @@ -146,6 +146,13 @@ def test_dwt_chisq_priors(): np.testing.assert_allclose(chisq, 1697.2230888243134) +def test_dwt_chisq_params_error(): + data = np.array([2.0, 0.0, 3.0, -2.0, -1.0, 2.0, 2.0, 0.0]) + model = np.ones(8) + params = np.array([1.0, 0.1]) + with pytest.raises(SystemExit): + chisq = ms.dwt_chisq(model, data, params) + def test_log_prior_uniform(): post = np.array([[3.0, 2.0], [3.1, 1.0], [3.6, 1.5]]) From 99d9a49297ef159a30433885c85884de4e5e322a Mon Sep 17 00:00:00 2001 From: Patricio Cubillos Date: Wed, 16 Oct 2019 16:46:32 +0200 Subject: [PATCH 3/6] Adjusted page numbering of plots.histogram() to match that of plot.trace(), i.e., start from zero. --- mc3/plots/plots.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mc3/plots/plots.py b/mc3/plots/plots.py index ed76e6b..266ece9 100644 --- a/mc3/plots/plots.py +++ b/mc3/plots/plots.py @@ -285,10 +285,10 @@ def histogram(posterior, pnames=None, thinning=1, fignum=300, ax.set_ylim(0, maxylim) if savefile is not None: - for j, fig in enumerate(figs): + for page, fig in enumerate(figs): if npages > 1: sf = os.path.splitext(savefile) - fig.savefig("{:s}_page{:02d}{:s}".format(sf[0], j+1, sf[1]), + fig.savefig("{:s}_page{:02d}{:s}".format(sf[0], page, sf[1]), bbox_inches='tight') else: fig.savefig(savefile, bbox_inches='tight') From cb63789684877236c7521812febb90202e747f6a Mon Sep 17 00:00:00 2001 From: Patricio Cubillos Date: Wed, 16 Oct 2019 16:51:37 +0200 Subject: [PATCH 4/6] Fixed bug with missing fignum number into plots.trace(). Set default fignums of figures to be in the 1000's. --- mc3/plots/plots.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/mc3/plots/plots.py b/mc3/plots/plots.py index 266ece9..fc03f0f 100644 --- a/mc3/plots/plots.py +++ b/mc3/plots/plots.py @@ -34,7 +34,7 @@ def trace(posterior, zchain=None, pnames=None, thinning=1, - burnin=0, fignum=100, savefile=None, fmt=".", ms=2.5, fs=11): + burnin=0, fignum=1000, savefile=None, fmt=".", ms=2.5, fs=11): """ Plot parameter trace MCMC sampling. @@ -98,7 +98,7 @@ def trace(posterior, zchain=None, pnames=None, thinning=1, axes = [] ipar = 0 for page in range(npages): - fig = plt.figure(page, figsize=(8.5,11.0)) + fig = plt.figure(fignum+page, figsize=(8.5,11.0)) plt.clf() plt.subplots_adjust(left=0.15, right=0.95, bottom=0.05, top=0.97, hspace=0.15) @@ -143,7 +143,7 @@ def trace(posterior, zchain=None, pnames=None, thinning=1, return axes -def histogram(posterior, pnames=None, thinning=1, fignum=300, +def histogram(posterior, pnames=None, thinning=1, fignum=1100, savefile=None, bestp=None, quantile=None, pdf=None, xpdf=None, ranges=None, axes=None, lw=2.0, fs=11, # Deprecated: Remove by 2020-07-01 @@ -296,7 +296,7 @@ def histogram(posterior, pnames=None, thinning=1, fignum=300, return axes -def pairwise(posterior, pnames=None, thinning=1, fignum=200, +def pairwise(posterior, pnames=None, thinning=1, fignum=1200, savefile=None, bestp=None, nbins=35, nlevels=20, absolute_dens=False, ranges=None, fs=11, rect=None, margin=0.01): """ @@ -453,7 +453,7 @@ def pairwise(posterior, pnames=None, thinning=1, fignum=200, def rms(binsz, rms, stderr, rmslo, rmshi, cadence=None, binstep=1, - timepoints=[], ratio=False, fignum=410, + timepoints=[], ratio=False, fignum=1300, yran=None, xran=None, savefile=None): """ Plot the RMS vs binsize curve. @@ -543,7 +543,7 @@ def rms(binsz, rms, stderr, rmslo, rmshi, cadence=None, binstep=1, def modelfit(data, uncert, indparams, model, nbins=75, - fignum=411, savefile=None, fmt="."): + fignum=1400, savefile=None, fmt="."): """ Plot the binned dataset with given uncertainties and model curves as a function of indparams. From dedea6e130d0f1aa41175762b941a8668cba7f3e Mon Sep 17 00:00:00 2001 From: Patricio Cubillos Date: Wed, 16 Oct 2019 18:16:52 +0200 Subject: [PATCH 5/6] Added extra condition to check best-fit vs MCMC lowest chisq, make sure that the best-fit values differ as well. --- mc3/sampler_driver.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/mc3/sampler_driver.py b/mc3/sampler_driver.py index 6f80f67..cf3206e 100644 --- a/mc3/sampler_driver.py +++ b/mc3/sampler_driver.py @@ -494,7 +494,8 @@ def sample(data=None, uncert=None, func=None, params=None, indparams=[], thinning, resume, log, **kwargs) if leastsq is not None: - if output['best_log_post']-fit_output['best_log_post'] > 3.0e-8: + if (output['best_log_post'] - fit_output['best_log_post'] > 5.0e-8 + and np.any((output['bestp'] - fit_output['bestp'])!=0.0)): log.warning("MCMC found a better fit than the minimizer:\n" "MCMC best-fitting parameters: (chisq={:.8g})\n{}\n" "Minimizer best-fitting parameters: (chisq={:.8g})\n{}". From 1877d37eba20dc337eaedd6f6a8445fba11d9ea4 Mon Sep 17 00:00:00 2001 From: Patricio Cubillos Date: Wed, 16 Oct 2019 18:24:49 +0200 Subject: [PATCH 6/6] Bumped MC3 version to 3.0.0 (no more pre-release, exclamation mark). Next merge fixes #116 --- mc3/VERSION.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mc3/VERSION.py b/mc3/VERSION.py index 4cb5e09..f5d4881 100644 --- a/mc3/VERSION.py +++ b/mc3/VERSION.py @@ -4,6 +4,6 @@ # MC3 Version: MC3_VER = 3 # Major version MC3_MIN = 0 # Minor version -MC3_REV = '0rc2' # Revision +MC3_REV = 0 # Revision __version__ = '{}.{}.{}'.format(MC3_VER, MC3_MIN, MC3_REV)