From 7854624837839576d9c82106fc27da1c84876686 Mon Sep 17 00:00:00 2001 From: springcoil Date: Wed, 21 Dec 2016 20:26:09 +0000 Subject: [PATCH 1/3] Fixing up tests Updating tests Updating tests again Added in one test for each sampler Updating tests Sampling errors Trying this commit Removing parameterized nose tests as they didn't work with the model context Moving asserts and removing arguments which are unnecessary - removing requirements.txt reference to nose_parameterized Changing shape Removing unnecessary test, and checking they run locally --- pymc3/tests/test_step.py | 50 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 49 insertions(+), 1 deletion(-) diff --git a/pymc3/tests/test_step.py b/pymc3/tests/test_step.py index b9b4ff9cb7..2aa343280e 100644 --- a/pymc3/tests/test_step.py +++ b/pymc3/tests/test_step.py @@ -3,16 +3,19 @@ from .checks import close_to from .models import simple_categorical, mv_simple, mv_simple_discrete, simple_2model +from .helpers import SeededTest +from pymc3 import df_summary, traceplot from pymc3.sampling import assign_step_methods, sample from pymc3.model import Model from pymc3.step_methods import (NUTS, BinaryGibbsMetropolis, CategoricalGibbsMetropolis, Metropolis, Slice, CompoundStep, MultivariateNormalProposal, HamiltonianMC) -from pymc3.distributions import Binomial, Normal, Bernoulli, Categorical +from pymc3.distributions import Binomial, Normal, Bernoulli, Categorical, Uniform from numpy.testing import assert_array_almost_equal import numpy as np from tqdm import tqdm +from scipy import stats class TestStepMethods(object): # yield test doesn't work subclassing unittest.TestCase @@ -238,3 +241,48 @@ def test_binomial(self): Binomial('x', 10, 0.5) steps = assign_step_methods(model, []) self.assertIsInstance(steps, Metropolis) + + +class TestSampleEstimates(SeededTest): + def test_posterior_estimate(self): + alpha_true, sigma_true = 1, 0.5 + beta_true = np.array([1, 2.5]) + + size = 100 + + X1 = np.random.randn(size) + X2 = np.random.randn(size) * 0.2 + Y = alpha_true + beta_true[0] * X1 + beta_true[1] * X2 + np.random.randn(size) * sigma_true + + + with Model() as model: + + alpha = Normal('alpha', mu=0, sd=10) + beta = Normal('beta', mu=0, sd=10, shape=2) + sigma = Uniform('sigma', lower=0.0, upper=1.0) + mu = alpha + beta[0] * X1 + beta[1] * X2 + Y_obs = Normal('Y_obs', mu=mu, sd=sigma, observed=Y) + + for step_method in (NUTS(), NUTS(is_cov=True), + Metropolis()): + trace = sample(100000, step=step_method, progressbar=False) + + # Using the Normal test from SciPy `scipy.stats.normaltest` + # See `https://github.com/scipy/scipy/blob/v0.18.1/scipy/stats/stats.py#L1352` + # `normaltest` returns a 2-tuple of the chi-squared statistic, and the associated + # p-value. Given the null hypothesis that `trace[-1000:]` came from the a normal + # distribution, the p-value represents the probability that a chi-squared statistic + # that large (or larger) would be seen. + + # If the p-value is very small, it means it's unlikely that the data came from a + # normal distribution. + # We define very small as 0.05, for the sake of this argument. + # We test this by evaluating if the p-value is above 0.05. + # We do the same for beta - using more burnin. + _, p_normal = stats.normaltest(trace[-3000:].alpha) + _, p_beta_0_normal = stats.normaltest(trace[-10000:].beta[:, 0]) + np.isclose(np.median(trace[-3000:].beta, 0), beta_true, rtol=0.1).all() + np.testing.assert_array_almost_equal(trace[-3000:].alpha, alpha_true, decimal=0) + np.testing.assert_array_almost_equal(np.median(trace[-3000:].sigma), sigma_true, decimal=1) + np.testing.assert_array_less(0.05, p_normal, verbose=True) + np.testing.assert_array_less(0.05, p_beta_0_normal, verbose=True) From c5c5a2d558b014e304e5491f3ef59f0466e67b57 Mon Sep 17 00:00:00 2001 From: Thomas Wiecki Date: Sat, 31 Dec 2016 17:59:44 +0100 Subject: [PATCH 2/3] TST Improve test. --- pymc3/tests/test_step.py | 64 ++++++++++++++++++---------------------- 1 file changed, 29 insertions(+), 35 deletions(-) diff --git a/pymc3/tests/test_step.py b/pymc3/tests/test_step.py index 2aa343280e..0c152b67c8 100644 --- a/pymc3/tests/test_step.py +++ b/pymc3/tests/test_step.py @@ -10,7 +10,7 @@ from pymc3.step_methods import (NUTS, BinaryGibbsMetropolis, CategoricalGibbsMetropolis, Metropolis, Slice, CompoundStep, MultivariateNormalProposal, HamiltonianMC) -from pymc3.distributions import Binomial, Normal, Bernoulli, Categorical, Uniform +from pymc3.distributions import Binomial, Normal, Bernoulli, Categorical, InverseGamma from numpy.testing import assert_array_almost_equal import numpy as np @@ -245,44 +245,38 @@ def test_binomial(self): class TestSampleEstimates(SeededTest): def test_posterior_estimate(self): - alpha_true, sigma_true = 1, 0.5 - beta_true = np.array([1, 2.5]) + alpha_true, sigma_true = 1., 0.5 + beta_true = 1. - size = 100 - - X1 = np.random.randn(size) - X2 = np.random.randn(size) * 0.2 - Y = alpha_true + beta_true[0] * X1 + beta_true[1] * X2 + np.random.randn(size) * sigma_true + size = 1000 + X = np.random.randn(size) + Y = alpha_true + beta_true * X + np.random.randn(size) * sigma_true + decimal = 1 with Model() as model: - - alpha = Normal('alpha', mu=0, sd=10) - beta = Normal('beta', mu=0, sd=10, shape=2) - sigma = Uniform('sigma', lower=0.0, upper=1.0) - mu = alpha + beta[0] * X1 + beta[1] * X2 + alpha = Normal('alpha', mu=0, sd=100, testval=alpha_true) + beta = Normal('beta', mu=0, sd=100, testval=beta_true) + sigma = InverseGamma('sigma', 10., testval=sigma_true) + mu = alpha + beta * X Y_obs = Normal('Y_obs', mu=mu, sd=sigma, observed=Y) - for step_method in (NUTS(), NUTS(is_cov=True), - Metropolis()): - trace = sample(100000, step=step_method, progressbar=False) - - # Using the Normal test from SciPy `scipy.stats.normaltest` - # See `https://github.com/scipy/scipy/blob/v0.18.1/scipy/stats/stats.py#L1352` - # `normaltest` returns a 2-tuple of the chi-squared statistic, and the associated - # p-value. Given the null hypothesis that `trace[-1000:]` came from the a normal - # distribution, the p-value represents the probability that a chi-squared statistic - # that large (or larger) would be seen. - - # If the p-value is very small, it means it's unlikely that the data came from a - # normal distribution. - # We define very small as 0.05, for the sake of this argument. - # We test this by evaluating if the p-value is above 0.05. + for step_method in (NUTS, Slice, Metropolis): + trace = sample(100000, step=step_method(), progressbar=False) + trace_ = trace[-3000:] + # We do the same for beta - using more burnin. - _, p_normal = stats.normaltest(trace[-3000:].alpha) - _, p_beta_0_normal = stats.normaltest(trace[-10000:].beta[:, 0]) - np.isclose(np.median(trace[-3000:].beta, 0), beta_true, rtol=0.1).all() - np.testing.assert_array_almost_equal(trace[-3000:].alpha, alpha_true, decimal=0) - np.testing.assert_array_almost_equal(np.median(trace[-3000:].sigma), sigma_true, decimal=1) - np.testing.assert_array_less(0.05, p_normal, verbose=True) - np.testing.assert_array_less(0.05, p_beta_0_normal, verbose=True) + np.testing.assert_almost_equal(np.mean(trace_.alpha), + alpha_true, decimal=decimal) + np.testing.assert_almost_equal(np.mean(trace_.beta), + beta_true, + decimal=decimal) + np.testing.assert_almost_equal(np.mean(trace_.sigma), + sigma_true, decimal=decimal) + + # Make sure posteriors are normal + _, p_alpha = stats.normaltest(trace_.alpha) + _, p_beta = stats.normaltest(trace_.beta) + # p-values should be > .05 to indiciate + np.testing.assert_array_less(0.05, p_alpha, verbose=True) + np.testing.assert_array_less(0.05, p_beta, verbose=True) From 043ac6bc6c4d147cd76ee42e3f3452e2e32827b0 Mon Sep 17 00:00:00 2001 From: Thomas Wiecki Date: Sat, 31 Dec 2016 19:18:14 +0100 Subject: [PATCH 3/3] TST Add thinning. --- pymc3/tests/test_step.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymc3/tests/test_step.py b/pymc3/tests/test_step.py index 0c152b67c8..7864068285 100644 --- a/pymc3/tests/test_step.py +++ b/pymc3/tests/test_step.py @@ -263,7 +263,7 @@ def test_posterior_estimate(self): for step_method in (NUTS, Slice, Metropolis): trace = sample(100000, step=step_method(), progressbar=False) - trace_ = trace[-3000:] + trace_ = trace[-300::5] # We do the same for beta - using more burnin. np.testing.assert_almost_equal(np.mean(trace_.alpha),