From ebaf7c476aedd74aca5f8fb11184947f2e441e73 Mon Sep 17 00:00:00 2001 From: Michael Osthege Date: Sun, 14 Nov 2021 20:10:23 +0100 Subject: [PATCH] Introduce continuous uni/multivariate model subclasses --- calibr8/__init__.py | 2 + calibr8/contrib/base.py | 2 +- calibr8/core.py | 243 ++++++----- calibr8/tests.py | 874 +++++++++++++++++++++------------------- 4 files changed, 576 insertions(+), 545 deletions(-) diff --git a/calibr8/__init__.py b/calibr8/__init__.py index ecd87df..586e4e1 100644 --- a/calibr8/__init__.py +++ b/calibr8/__init__.py @@ -14,6 +14,8 @@ from .core import ( CalibrationModel, DistributionMixin, + ContinuousMultivariateModel, + ContinuousUnivariateModel, InferenceResult, ContinuousMultivariateInference, ContinuousUnivariateInference, diff --git a/calibr8/contrib/base.py b/calibr8/contrib/base.py index da8518f..5c55874 100644 --- a/calibr8/contrib/base.py +++ b/calibr8/contrib/base.py @@ -8,7 +8,7 @@ from .. import core -class BaseModelT(core.CalibrationModel, noise.StudentTNoise): +class BaseModelT(core.ContinuousUnivariateModel, noise.StudentTNoise): pass diff --git a/calibr8/core.py b/calibr8/core.py index 517f042..76c8b5e 100644 --- a/calibr8/core.py +++ b/calibr8/core.py @@ -272,101 +272,6 @@ def hdi_objective(x): return hdi_lower, hdi_upper -def _infer_univariate_independent( - likelihood, - y:Union[int, float, numpy.ndarray], - *, - lower:float, - upper:float, - steps:int=300, - ci_prob:float=1, -) -> ContinuousUnivariateInference: - """Infer the posterior distribution of the independent variable given the observations of the dependent variable. - The calculation is done numerically by integrating the likelihood in a certain interval [upper,lower]. - This is identical to the posterior with a Uniform (lower,upper) prior. If precentiles are provided, the interval of - the PDF will be shortened. - - Parameters - ---------- - likelihood : callable - A likelihood function to integrate over. - Needs to accept keyword arguments: x, y, scan_x - y : int, float, array - one or more observations at the same x - lower : float - lower limit for uniform distribution of prior - upper : float - upper limit for uniform distribution of prior - steps : int - steps between lower and upper or steps between the percentiles (default 300) - ci_prob : float - Probability level for ETI and HDI credible intervals. - If 1 (default), the complete interval [upper,lower] will be returned, - else the PDFs will be trimmed to the according probability interval; - float must be in the interval (0,1] - - Returns - ------- - posterior : ContinuousUnivariateInference - the result of the numeric posterior calculation - """ - y = numpy.atleast_1d(y) - - likelihood_integral, _ = scipy.integrate.quad( - func=lambda x: likelihood(x=x, y=y), - # by restricting the integral into the interval [a,b], the resulting PDF is - # identical to the posterior with a Uniform(a, b) prior. - # 1. prior probability is constant in [a,b] - # 2. prior probability is 0 outside of [a,b] - # > numerical integral is only computed in [a,b], but because of 1. and 2., it's - # identical to the integral over [-∞,+∞] - a=lower, b=upper, - ) - - # high resolution x-coordinates for integration - # the first integration is just to find the peak - x_integrate = numpy.linspace(lower, upper, 10_000) - area = scipy.integrate.cumtrapz(likelihood(x=x_integrate, y=y, scan_x=True), x_integrate, initial=0) - cdf = area / area[-1] - - # now we find a high-resolution CDF for (1-shrink) of the probability mass - shrink = 0.00001 - xfrom, xto = _get_eti(x_integrate, cdf, 1 - shrink) - x_integrate = numpy.linspace(xfrom, xto, 100_000) - area = scipy.integrate.cumtrapz(likelihood(x=x_integrate, y=y, scan_x=True), x_integrate, initial=0) - cdf = (area / area[-1]) * (1 - shrink) + shrink / 2 - - # TODO: create a smart x-vector from the CDF with varying stepsize - - if ci_prob != 1: - if not isinstance(ci_prob, (int, float)) or not (0 < ci_prob <= 1): - raise ValueError(f'Unexpected `ci_prob` value of {ci_prob}. Expected float in interval (0, 1].') - - # determine the interval bounds from the high-resolution CDF - eti_lower, eti_upper = _get_eti(x_integrate, cdf, ci_prob) - hdi_lower, hdi_upper = _get_hdi(x_integrate, cdf, ci_prob, eti_lower, eti_upper, history=None) - - eti_x = numpy.linspace(eti_lower, eti_upper, steps) - hdi_x = numpy.linspace(hdi_lower, hdi_upper, steps) - eti_pdf = likelihood(x=eti_x, y=y, scan_x=True) / likelihood_integral - hdi_pdf = likelihood(x=hdi_x, y=y, scan_x=True) / likelihood_integral - eti_prob = _interval_prob(x_integrate, cdf, eti_lower, eti_upper) - hdi_prob = _interval_prob(x_integrate, cdf, hdi_lower, hdi_upper) - else: - x = numpy.linspace(lower, upper, steps) - eti_x = hdi_x = x - eti_pdf = hdi_pdf = likelihood(x=x, y=y, scan_x=True) / likelihood_integral - eti_prob = hdi_prob = 1 - - median = x_integrate[numpy.argmin(numpy.abs(cdf - 0.5))] - - return ContinuousUnivariateInference( - median, - eti_x=eti_x, eti_pdf=eti_pdf, eti_prob=eti_prob, - hdi_x=hdi_x, hdi_pdf=hdi_pdf, hdi_prob=hdi_prob, - ) - - class DistributionMixin: """Maps the values returned by `CalibrationModel.predict_dependent` to a SciPy distribution and its parameters, and optionally also to @@ -397,7 +302,14 @@ def _inherits_noisemodel(cls): class CalibrationModel(DistributionMixin): """A parent class providing the general structure of a calibration model.""" - def __init__(self, independent_key:str, dependent_key:str, *, theta_names:typing.Tuple[str], ndim=1): + def __init__( + self, + independent_key:str, + dependent_key:str, + *, + theta_names:typing.Tuple[str], + ndim: int, + ): """Creates a CalibrationModel object. Parameters @@ -409,9 +321,7 @@ def __init__(self, independent_key:str, dependent_key:str, *, theta_names:typing theta_names : optional, tuple of str names of the model parameters ndim : int - Number of dimensions that the independent variable has. - Most calibrations are univariate (ndim=1). - Multivariate calibrations (ndim > 1) may need to override the `infer_independent` implementation. + Number of independent dimensions in the model. """ if not _inherits_noisemodel(type(self)): warnings.warn( @@ -512,13 +422,11 @@ def infer_independent( self, y:Union[int,float,numpy.ndarray], *, - lower:float, upper:float, steps:int=300, - ci_prob:float=1 + lower, + upper, ) -> InferenceResult: """Infer the posterior distribution of the independent variable given the observations of the dependent variable. - The calculation is done numerically by integrating the likelihood in a certain interval [upper,lower]. - This is identical to the posterior with a Uniform (lower,upper) prior. If precentiles are provided, the interval of - the PDF will be shortened. + A Uniform (lower,upper) prior is applied. Parameters ---------- @@ -528,33 +436,13 @@ def infer_independent( lower limit for uniform distribution of prior upper : float upper limit for uniform distribution of prior - steps : int - steps between lower and upper or steps between the percentiles (default 300) - ci_prob : float - Probability level for ETI and HDI credible intervals. - If 1 (default), the complete interval [upper,lower] will be returned, - else the PDFs will be trimmed to the according probability interval; - float must be in the interval (0,1] Returns ------- posterior : InferenceResult - the result of the numeric posterior calculation + Result of the independent variable inference. """ - if self.ndim != 1: - raise NotImplementedError( - f"This calibration model seems to be multivariate (ndim={self.ndim})," - " but does not have an override implementation of .infer_independent()." - ) - y = numpy.atleast_1d(y) - return _infer_univariate_independent( - self.likelihood, - y, - lower=lower, - upper=upper, - steps=steps, - ci_prob=ci_prob, - ) + raise NotImplementedError(f"This calibration model does not implement an .infer_independent() method.") def loglikelihood( self, @@ -786,6 +674,109 @@ def load(cls, filepath: os.PathLike): return obj +class ContinuousUnivariateModel(CalibrationModel): + def __init__(self, independent_key: str, dependent_key: str, *, theta_names: typing.Tuple[str]): + super().__init__(independent_key, dependent_key, theta_names=theta_names, ndim=1) + + def infer_independent( + self, + y:Union[int, float, numpy.ndarray], + *, + lower:float, + upper:float, + steps:int=300, + ci_prob:float=1, + ) -> ContinuousUnivariateInference: + """Infer the posterior distribution of the independent variable given the observations of the dependent variable. + The calculation is done numerically by integrating the likelihood in a certain interval [upper,lower]. + This is identical to the posterior with a Uniform (lower,upper) prior. + + Parameters + ---------- + y : int, float, array + One or more observations at the same x. + lower : float + Lower limit for uniform distribution of prior. + upper : float + Upper limit for uniform distribution of prior. + steps : int + Steps between lower and upper or steps between the percentiles (default 300). + ci_prob : float + Probability level for ETI and HDI credible intervals. + If 1 (default), the complete interval [upper,lower] will be returned, + else the PDFs will be trimmed to the according probability interval; + float must be in the interval (0,1] + + Returns + ------- + posterior : ContinuousUnivariateInference + Result of the numeric posterior calculation. + """ + y = numpy.atleast_1d(y) + + likelihood_integral, _ = scipy.integrate.quad( + func=lambda x: self.likelihood(x=x, y=y), + # by restricting the integral into the interval [a,b], the resulting PDF is + # identical to the posterior with a Uniform(a, b) prior. + # 1. prior probability is constant in [a,b] + # 2. prior probability is 0 outside of [a,b] + # > numerical integral is only computed in [a,b], but because of 1. and 2., it's + # identical to the integral over [-∞,+∞] + a=lower, b=upper, + ) + + # high resolution x-coordinates for integration + # the first integration is just to find the peak + x_integrate = numpy.linspace(lower, upper, 10_000) + area = scipy.integrate.cumtrapz(self.likelihood(x=x_integrate, y=y, scan_x=True), x_integrate, initial=0) + cdf = area / area[-1] + + # now we find a high-resolution CDF for (1-shrink) of the probability mass + shrink = 0.00001 + xfrom, xto = _get_eti(x_integrate, cdf, 1 - shrink) + x_integrate = numpy.linspace(xfrom, xto, 100_000) + area = scipy.integrate.cumtrapz(self.likelihood(x=x_integrate, y=y, scan_x=True), x_integrate, initial=0) + cdf = (area / area[-1]) * (1 - shrink) + shrink / 2 + + # TODO: create a smart x-vector from the CDF with varying stepsize + + if ci_prob != 1: + if not isinstance(ci_prob, (int, float)) or not (0 < ci_prob <= 1): + raise ValueError(f'Unexpected `ci_prob` value of {ci_prob}. Expected float in interval (0, 1].') + + # determine the interval bounds from the high-resolution CDF + eti_lower, eti_upper = _get_eti(x_integrate, cdf, ci_prob) + hdi_lower, hdi_upper = _get_hdi(x_integrate, cdf, ci_prob, eti_lower, eti_upper, history=None) + + eti_x = numpy.linspace(eti_lower, eti_upper, steps) + hdi_x = numpy.linspace(hdi_lower, hdi_upper, steps) + eti_pdf = self.likelihood(x=eti_x, y=y, scan_x=True) / likelihood_integral + hdi_pdf = self.likelihood(x=hdi_x, y=y, scan_x=True) / likelihood_integral + eti_prob = _interval_prob(x_integrate, cdf, eti_lower, eti_upper) + hdi_prob = _interval_prob(x_integrate, cdf, hdi_lower, hdi_upper) + else: + x = numpy.linspace(lower, upper, steps) + eti_x = hdi_x = x + eti_pdf = hdi_pdf = self.likelihood(x=x, y=y, scan_x=True) / likelihood_integral + eti_prob = hdi_prob = 1 + + median = x_integrate[numpy.argmin(numpy.abs(cdf - 0.5))] + + return ContinuousUnivariateInference( + median, + eti_x=eti_x, eti_pdf=eti_pdf, eti_prob=eti_prob, + hdi_x=hdi_x, hdi_pdf=hdi_pdf, hdi_prob=hdi_prob, + ) + + +class ContinuousMultivariateModel(CalibrationModel): + def __init__(self, independent_key: str, dependent_key: str, *, theta_names: typing.Tuple[str], ndim: int): + super().__init__(independent_key, dependent_key, theta_names=theta_names, ndim=ndim) + + def infer_independent(self, y: Union[int, float, numpy.ndarray], *, lower, upper) -> ContinuousMultivariateInference: + return super().infer_independent(y, lower=lower, upper=upper) + + def logistic(x, theta): """4-parameter logistic model. diff --git a/calibr8/tests.py b/calibr8/tests.py index 0fb073d..3aab7f4 100644 --- a/calibr8/tests.py +++ b/calibr8/tests.py @@ -41,10 +41,10 @@ class _TestModel(calibr8.CalibrationModel, calibr8.NormalNoise): def __init__(self, independent_key=None, dependent_key=None, theta_names=None): if theta_names is None: theta_names = tuple('a,b,c'.split(',')) - super().__init__(independent_key='I', dependent_key='D', theta_names=theta_names) + super().__init__(independent_key='I', dependent_key='D', theta_names=theta_names, ndim=1) -class _TestBivariateLinearModel(calibr8.CalibrationModel, calibr8.NormalNoise): +class _TestBivariateLinearModel(calibr8.ContinuousMultivariateModel, calibr8.NormalNoise): def __init__(self, independent_key: str=None, dependent_key: str=None, theta_names=None): super().__init__(independent_key="x1,x2", dependent_key="y", theta_names="i,s1,s2,sd".split(","), ndim=2) @@ -201,18 +201,18 @@ def __init__(self): pass def test_constructor_signature_check(self): - class EM_OK(calibr8.CalibrationModel): + class EM_OK(calibr8.ContinuousUnivariateModel): def __init__(self, arg1=1, *, kwonly=2, kwonlydefault=4): super().__init__('I', 'D', theta_names=tuple('abc')) EM_OK() - class EM_args(calibr8.CalibrationModel): + class EM_args(calibr8.ContinuousUnivariateModel): def __init__(self, arg1): super().__init__('I', 'D', theta_names=tuple('abc')) with pytest.raises(TypeError, match=r"constructor must not have any required \(kw\)arguments"): EM_args(arg1=3) - class EM_kwargs(calibr8.CalibrationModel): + class EM_kwargs(calibr8.ContinuousUnivariateModel): def __init__(self, *, kwonly, kwonlydefault=4): super().__init__('I', 'D', theta_names=tuple('abc')) with pytest.raises(TypeError, match=r"constructor must not have any required \(kw\)arguments"): @@ -230,11 +230,7 @@ def test_exceptions(self): cmodel.predict_independent(x) with pytest.raises(NotImplementedError, match="predict_dependent function"): cmodel.loglikelihood(y=y, x=x, theta=[1, 2, 3]) - with pytest.raises(ValueError, match=r"Unexpected `ci_prob`"): - cmodel.loglikelihood = lambda x, y, theta: 1 - cmodel.infer_independent(y, lower=0, upper=10, steps=10, ci_prob=None) - with pytest.raises(NotImplementedError, match="seems to be multivariate"): - cmodel.ndim = 3 + with pytest.raises(NotImplementedError, match=r"does not implement an .infer_independent\(\)"): cmodel.infer_independent(2, lower=0, upper=5) pass @@ -362,6 +358,248 @@ def test_likelihood_multivariate(self): pass +class TestContinuousUnivariateModel: + def test_continuous_univariate_exceptions(self): + cmodel = _TestLogisticModel(scale_degree=1) + cmodel.theta_fitted = [0, 4, 2, 1, 1, 0, 2, 1.4] + with pytest.raises(ValueError, match=r"Unexpected `ci_prob`"): + cmodel.infer_independent(42, lower=0, upper=10, steps=10, ci_prob=None) + pass + + @pytest.mark.skipif(not HAS_PYMC, reason='requires PyMC') + def test_symbolic_loglikelihood_checks_and_warnings(self): + cmodel = _TestPolynomialModel(independent_key='S', dependent_key='A', mu_degree=1, scale_degree=1) + cmodel.theta_fitted = [0, 1, 0.1, 1, 5] + + # create test data + x_true = numpy.array([1,2,3,4,5]) + y_obs = cmodel.predict_dependent(x_true)[0] + + with pm.Model(): + x_hat = pm.Uniform("x_hat", shape=5) + with pytest.raises(ValueError, match="`name` must be specified"): + cmodel.loglikelihood(x=x_hat, y=y_obs) + + with pytest.warns(DeprecationWarning, match="Use `name` instead"): + cmodel.loglikelihood(x=x_hat, y=y_obs, replicate_id='A01', dependent_key='A') + pass + + @pytest.mark.skipif(not HAS_PYMC, reason='requires PyMC') + def test_symbolic_loglikelihood(self): + cmodel = _TestPolynomialModel(independent_key='S', dependent_key='A', mu_degree=1, scale_degree=1) + cmodel.theta_fitted = [0, 1, 0.1, 1, 3] + + # create test data + x_true = numpy.array([1,2,3,4,5]) + y_obs = cmodel.predict_dependent(x_true)[0] + + x_hat = at.vector() + x_hat.tag.test_value = x_true + L = cmodel.loglikelihood(x=x_hat, y=y_obs, name='L_A01_A') + assert isinstance(L, at.TensorVariable) + assert L.ndim == 0 + + # compare the two loglikelihood computation methods + x_test = numpy.random.normal(x_true, scale=0.1) + actual = L.eval({ + x_hat: x_test + }) + expected = cmodel.loglikelihood(x=x_test, y=y_obs) + assert numpy.ndim(expected) == 0 + assert numpy.ndim(actual) == 0 + numpy.testing.assert_almost_equal(actual, expected, 6) + pass + + @pytest.mark.skipif(not HAS_PYMC, reason='requires PyMC') + def test_symbolic_loglikelihood_in_modelcontext(self): + cmodel = _TestPolynomialModel(independent_key='S', dependent_key='A', mu_degree=1, scale_degree=1) + cmodel.theta_fitted = [0, 0.5, 0.1, 1, 7] + + # create test data + x_true = numpy.array([1,2,3,4,5]) + y_obs = cmodel.predict_dependent(x_true)[0] + + # create a PyMC model using the calibration model + with pm.Model() as pmodel: + x_hat = pm.Uniform("x_hat", 0, 1, shape=x_true.shape, transform=None) + L = cmodel.loglikelihood(x=x_hat, y=y_obs, name='L_A01_A') + assert isinstance(L, at.TensorVariable) + assert L.ndim == 0 + + # PyMC v4 returns the RV, but for .eval() we need the RV-value-variable + if pm.__version__[0] != "3": + x_hat = pmodel.rvs_to_values[x_hat] + + # compare the two loglikelihood computation methods + x_test = numpy.random.normal(x_true, scale=0.1) + actual = L.eval({ + x_hat: x_test + }) + expected = cmodel.loglikelihood(x=x_test, y=y_obs) + assert numpy.ndim(expected) == 0 + assert numpy.ndim(actual) == 0 + numpy.testing.assert_almost_equal(actual, expected, 6) + pass + + @pytest.mark.parametrize("x", [ + numpy.array([1,2,3]), + 4, + ]) + @pytest.mark.parametrize("y", [ + numpy.array([2,4,8]), + 5, + ]) + def test_loglikelihood(self, x, y): + cmodel = _TestPolynomialModel(independent_key='S', dependent_key='OD', mu_degree=1, scale_degree=1) + cmodel.theta_fitted = [0, 1, 0.1, 1.6, 3] + + actual = cmodel.loglikelihood(y=y, x=x) + assert numpy.ndim(actual) == 0 + mu, scale, df = cmodel.predict_dependent(x, theta=cmodel.theta_fitted) + expected = numpy.sum(stats.t.logpdf(x=y, loc=mu, scale=scale, df=df)) + numpy.testing.assert_equal(actual, expected) + return + + def test_loglikelihood_exceptions(self): + cmodel = _TestPolynomialModel(independent_key='S', dependent_key='OD', mu_degree=1, scale_degree=1) + with pytest.raises(Exception, match="No parameter vector"): + cmodel.loglikelihood(y=[2,3], x=[4,5]) + + cmodel.theta_fitted = [0, 1, 0.1, 1.6, 2] + + with pytest.raises(TypeError): + cmodel.loglikelihood(4, x=[2,3]) + with pytest.raises(ValueError, match="Input x must be"): + cmodel.loglikelihood(y=[2,3], x="hello") + with pytest.raises(ValueError, match="Input y must be"): + cmodel.loglikelihood(y="🤔", x=2) + with pytest.raises(ValueError, match="operands could not be broadcast"): + cmodel.loglikelihood(y=[1,2,3], x=[1,2]) + return + + def test_likelihood(self): + # use a linear model with intercept 1 and slope 0.5 + cmodel = _TestPolynomialModel(independent_key="I", dependent_key="D", mu_degree=1) + cmodel.theta_fitted = [1, 0.5, 0.5, 4] + + assert numpy.isscalar(cmodel.likelihood(y=2, x=3)) + assert numpy.isscalar(cmodel.likelihood(y=[2, 3], x=[3, 4])) + + with pytest.raises(ValueError, match="operands could not be broadcast"): + cmodel.likelihood(y=[1,2,3], x=[1,2]) + + x_dense = numpy.linspace(0, 4, 501) + actual = cmodel.likelihood(x=x_dense, y=2, scan_x=True) + assert numpy.ndim(actual) == 1 + # the maximum likelihood should be at x=2 + assert x_dense[numpy.argmax(actual)] == 2 + pass + + @pytest.mark.xfail(reason="Draft test case, see https://github.com/JuBiotech/calibr8/issues/15") + def test_likelihood_nobroadcasting_fallback(self): + class _TestSwitchableBroadcastingModel(calibr8.CalibrationModel, calibr8.NormalNoise): + def __init__(self): + self.x_shapes = [] + self.can_broadcast = False + super().__init__(independent_key="I", dependent_key="D", theta_names="i,s,sd".split(","), ndim=1) + + def predict_dependent(self, x, *, theta=None): + if theta is None: + theta = self.theta_fitted + i, s, sd = self.theta_fitted + # This part broadcasts just fine + x = numpy.array(x) + mu = i + x * s + return mu, sd + + def loglikelihood(self, *, y, x, **kwargs): + # This overrides the native CalibrationModel.loglikelihood with one + # that can be externally set to return non-broadcasted results. + LL_broadcasted = super().loglikelihood(y=y, x=x, **kwargs) + if not self.can_broadcast: + return LL_broadcasted.sum() + return LL_broadcasted + + cm = _TestSwitchableBroadcastingModel() + cm.theta_fitted = [0.5, 0.6, 0.7] + + X = numpy.array([1, 2, 3]) + Y = [1.1, 1.7, 2.3] + y_obs = [1.0, 1.7, 2.4] + + # Check the prediction of the test model + for x, y in zip(X, Y): + assert cm.predict_dependent(x) == (y, 0.7) + + # The predict_dependent can broadcast + mu, sd = cm.predict_dependent(X) + numpy.testing.assert_array_equal(mu, Y) + assert sd == 0.7 + + # Test the switching between can_broadcast modes + cm.can_broadcast = True + LL = cm.loglikelihood(x=X[..., None], y=y_obs) + assert numpy.shape(LL) == (3,) + + cm.can_broadcast = False + LL = cm.loglikelihood(x=X[..., None], y=y_obs) + assert numpy.shape(LL) == () + + # The CalibrationModel.likelihood should give the same results either way. + cm.can_broadcast = True + L_broadcasted = cm.likelihood(x=X[..., None], y=y_obs, scan_x=True) + cm.can_broadcast = False + L_looped = cm.likelihood(x=X[..., None], y=y_obs, scan_x=True) + + numpy.testing.assert_array_equal(L_broadcasted, L_looped) + + # Of course the values should also be correct. + L_expected = numpy.exp([ + cm.scipy_dist.logpdf(x=y_obs, loc=mui, scale=sd).sum() + for mui in mu + ]) + assert numpy.shape(L_expected) == (3,) + numpy.testing.assert_array_equal(L_broadcasted, L_expected) + pass + + def test_infer_independent(self): + em = _TestPolynomialModel(independent_key='S', dependent_key='A365', mu_degree=1, scale_degree=1) + em.theta_fitted = [0, 2, 0.1, 1, 3] + pst = em.infer_independent(y=1, lower=0, upper=20, steps=876) + + assert len(pst.eti_x) == len(pst.eti_pdf) + assert len(pst.hdi_x) == len(pst.hdi_pdf) + assert tuple(pst.eti_x[[0, -1]]) == (0, 20) + assert tuple(pst.hdi_x[[0, -1]]) == (0, 20) + assert (numpy.isclose(scipy.integrate.cumtrapz(pst.eti_pdf, pst.eti_x)[-1], 1, atol=0.0001)) + assert (numpy.isclose(scipy.integrate.cumtrapz(pst.hdi_pdf, pst.hdi_x)[-1], 1, atol=0.0001)) + assert pst.eti_lower == pst.hdi_lower == 0 + assert pst.eti_upper == pst.hdi_upper == 20 + assert pst.eti_prob == 1 + assert pst.hdi_prob == 1 + + # check trimming to [2.5,97.5] interval + pst = em.infer_independent(y=[1, 2], lower=0, upper=20, steps=1775, ci_prob=0.95) + + assert len(pst.eti_x) == len(pst.eti_pdf) + assert len(pst.hdi_x) == len(pst.hdi_pdf) + assert numpy.isclose(pst.eti_prob, 0.95, atol=0.0001) + assert numpy.isclose(pst.hdi_prob, 0.95, atol=0.0001) + assert numpy.isclose(scipy.integrate.cumtrapz(pst.eti_pdf, pst.eti_x)[-1], 0.95, atol=0.0001) + assert numpy.isclose(scipy.integrate.cumtrapz(pst.hdi_pdf, pst.hdi_x)[-1], 0.95, atol=0.0001) + assert pst.eti_lower == pst.eti_x[0] + assert pst.eti_upper == pst.eti_x[-1] + assert pst.hdi_lower == pst.hdi_x[0] + assert pst.hdi_upper == pst.hdi_x[-1] + + # check that error are raised by wrong input + with pytest.raises(ValueError): + em.infer_independent(y=1, lower=0, upper=20, steps=1000, ci_prob=(-1)) + with pytest.raises(ValueError): + em.infer_independent(y=1, lower=0, upper=20, steps=1000, ci_prob=(97.5)) + pass + + class TestUnivariateInferenceHelpers: def test_get_eti(self): mu = 0.2 @@ -432,26 +670,72 @@ def test_get_hdi(self, initial_guess): pass -class TestModelFunctions: - def test_logistic(self): - x = numpy.array([1.,2.,4.]) - theta = [2,2,4,1] - expected = 2*2-4+(2*(4-2))/(1+numpy.exp(-2*1/(4-2)*(x-2))) - true = calibr8.logistic(x, theta) - assert (numpy.array_equal(true, expected)) - return - - def test_inverse_logistic(self): - x = numpy.array([1.,2.,4.]) - theta = [2,2,4,1] - forward = calibr8.logistic(x, theta) - reverse = calibr8.inverse_logistic(forward, theta) - assert (numpy.allclose(x, reverse)) - return +class TestContinuousMultivariateModel: + def test_likelihood_multivariate(self): + cm = _TestBivariateLinearModel() + cm.theta_fitted = (0.5, 1, 2, 0.4) - def test_asymmetric_logistic(self): - L_L = -3 - L_U = 4 + X = [ + [1, 2], + [0.5, 1], + [1.5, 2.5], + ] + Y = [5.5, 3.0, 7.0] + y_obs = [5.6, 3.0, 6.9] + + # Check the expected output for one coordinate + assert cm.predict_dependent(X[0]) == (Y[0], 0.4) + + # And its broadcasting for multiple coordinates + numpy.testing.assert_array_equal(cm.predict_dependent(X)[0], Y) + + # Expected likelihoods of coordinate/observation paris + LL_elemwise = cm.scipy_dist.logpdf(x=y_obs, loc=Y, scale=0.4).sum() + + # Expected likelihoods of all observation at each coordinate + LL_scan = numpy.array([ + cm.scipy_dist.logpdf(x=y_obs, loc=y, scale=0.4).sum() + for y in Y + ]) + + assert LL_elemwise.shape == () + assert LL_scan.shape == (3,) + + # Testing the underlying loglikelihood + numpy.testing.assert_array_equal(cm.loglikelihood(y=y_obs, x=X), LL_elemwise) + + # Now check if likelihood wraps it correctly + numpy.testing.assert_array_equal(cm.likelihood(y=y_obs, x=X, scan_x=False), numpy.exp(LL_elemwise)) + numpy.testing.assert_array_equal(cm.likelihood(y=y_obs, x=X, scan_x=True), numpy.exp(LL_scan)) + pass + + def test_infer_independent_not_implemented(self): + cm = _TestBivariateLinearModel() + with pytest.raises(NotImplementedError, match=r"does not implement an .infer_independent\(\)"): + cm.infer_independent(y=1, lower=[0, 0], upper=[1, 1]) + pass + + +class TestModelFunctions: + def test_logistic(self): + x = numpy.array([1.,2.,4.]) + theta = [2,2,4,1] + expected = 2*2-4+(2*(4-2))/(1+numpy.exp(-2*1/(4-2)*(x-2))) + true = calibr8.logistic(x, theta) + assert (numpy.array_equal(true, expected)) + return + + def test_inverse_logistic(self): + x = numpy.array([1.,2.,4.]) + theta = [2,2,4,1] + forward = calibr8.logistic(x, theta) + reverse = calibr8.inverse_logistic(forward, theta) + assert (numpy.allclose(x, reverse)) + return + + def test_asymmetric_logistic(self): + L_L = -3 + L_U = 4 I_x = 4.5 S = 3.3 c = -1 @@ -782,395 +1066,6 @@ def test_plot_model_band_xlim(self): pass -class TestContribBase: - def test_cant_instantiate_base_models(self): - with pytest.raises(TypeError): - calibr8.BaseModelT(independent_key='I', dependent_key='D') - with pytest.raises(TypeError): - calibr8.BaseAsymmetricLogisticT(independent_key='I', dependent_key='D') - with pytest.raises(TypeError): - calibr8.BasePolynomialModelT(independent_key='I', dependent_key='D', mu_degree=1, scale_degree=1) - pass - - def test_base_polynomial_t(self): - for mu_degree, scale_degree in [(1, 0), (1, 1), (1, 2), (2, 0)]: - theta_mu = (2.2, 1.2, 0.2)[:mu_degree+1] - theta_scale = (3.1, 0.4, 0.2)[:scale_degree+1] - theta = theta_mu + theta_scale + (1,) - - em = _TestPolynomialModel(independent_key='I', dependent_key='D', mu_degree=mu_degree, scale_degree=scale_degree) - assert len(em.theta_names) == mu_degree+1 + scale_degree+1 + 1 - assert len(em.theta_names) == len(theta) - - x = numpy.linspace(0, 10, 3) - mu, scale, df = em.predict_dependent(x, theta=theta) - - - expected = numpy.polyval(theta_mu[::-1], x) - numpy.testing.assert_array_equal(mu, expected) - - expected = numpy.polyval(theta_scale[::-1], mu) - numpy.testing.assert_array_equal(scale, expected) - - numpy.testing.assert_array_equal(df, 1) - pass - - def test_base_polynomial_t_inverse(self): - for mu_degree, scale_degree in [(1, 0), (1, 1), (1, 2), (2, 0)]: - theta_mu = (2.2, 1.2, 0.2)[:mu_degree+1] - theta_scale = (3.1, 0.4, 0.2)[:scale_degree+1] - theta = theta_mu + theta_scale + (1,) - - em = _TestPolynomialModel(independent_key='I', dependent_key='D', mu_degree=mu_degree, scale_degree=scale_degree) - em.theta_fitted = theta - assert len(em.theta_names) == mu_degree+1 + scale_degree+1 + 1 - assert len(em.theta_names) == len(theta) - - x = numpy.linspace(0, 10, 7) - mu, scale, df = em.predict_dependent(x, theta=theta) - - if mu_degree < 2: - x_inverse = em.predict_independent(mu) - numpy.testing.assert_array_almost_equal(x_inverse, x) - else: - with pytest.raises(NotImplementedError): - em.predict_independent(mu) - pass - - def test_base_asymmetric_logistic_t(self): - for scale_degree in [0, 1, 2]: - theta_mu = (-0.5, 0.5, 1, 1, -1) - theta_scale = (3.1, 0.4, 0.2)[:scale_degree+1] - theta = theta_mu + theta_scale + (1,) - - em = _TestLogisticModel(independent_key='I', dependent_key='D', scale_degree=scale_degree) - assert len(em.theta_names) == 5 + scale_degree+1 + 1 - assert len(em.theta_names) == len(theta) - - x = numpy.linspace(1, 5, 3) - mu, scale, df = em.predict_dependent(x, theta=theta) - - expected = calibr8.asymmetric_logistic(x, theta_mu) - numpy.testing.assert_array_equal(mu, expected) - - expected = numpy.polyval(theta_scale[::-1], mu) - numpy.testing.assert_array_equal(scale, expected) - - numpy.testing.assert_array_equal(df, 1) - pass - - def test_base_asymmetric_logistic_t_inverse(self): - for scale_degree in [0, 1, 2]: - theta_mu = (-0.5, 0.5, 1, 1, -1) - theta_scale = (3.1, 0.4, 0.2)[:scale_degree+1] - theta = theta_mu + theta_scale + (1,) - - em = _TestLogisticModel(independent_key='I', dependent_key='D', scale_degree=scale_degree) - em.theta_fitted = theta - - x = numpy.linspace(1, 5, 7) - mu, scale, df = em.predict_dependent(x, theta=theta) - - x_inverse = em.predict_independent(mu) - numpy.testing.assert_array_almost_equal(x_inverse, x) - pass - - def test_base_xlog_asymmetric_logistic_t(self): - for scale_degree in [0, 1, 2]: - theta_mu = (-0.5, 0.5, 1, 1, -1) - theta_scale = (3.1, 0.4, 0.2)[:scale_degree+1] - theta = theta_mu + theta_scale + (1,) - - em = _TestLogIndependentLogisticModel(independent_key='I', dependent_key='D', scale_degree=scale_degree) - assert len(em.theta_names) == 5 + scale_degree+1 + 1 - assert len(em.theta_names) == len(theta) - - x = numpy.linspace(1, 5, 3) - mu, scale, df = em.predict_dependent(x, theta=theta) - - expected = calibr8.xlog_asymmetric_logistic(x, theta_mu) - numpy.testing.assert_array_equal(mu, expected) - - expected = numpy.polyval(theta_scale[::-1], mu) - numpy.testing.assert_array_equal(scale, expected) - - numpy.testing.assert_array_equal(df, 1) - pass - - def test_base_xlog_asymmetric_logistic_t_inverse(self): - for scale_degree in [0, 1, 2]: - theta_mu = (-0.5, 0.5, 1, 1, -1) - theta_scale = (3.1, 0.4, 0.2)[:scale_degree+1] - theta = theta_mu + theta_scale + (1,) - - em = _TestLogIndependentLogisticModel(independent_key='I', dependent_key='D', scale_degree=scale_degree) - em.theta_fitted = theta - - x = numpy.linspace(1, 5, 7) - mu, scale, df = em.predict_dependent(x, theta=theta) - - x_inverse = em.predict_independent(mu) - numpy.testing.assert_array_almost_equal(x_inverse, x) - pass - - -class TestBasePolynomialModelT: - def test_infer_independent(self): - em = _TestPolynomialModel(independent_key='S', dependent_key='A365', mu_degree=1, scale_degree=1) - em.theta_fitted = [0, 2, 0.1, 1, 3] - pst = em.infer_independent(y=1, lower=0, upper=20, steps=876) - - assert len(pst.eti_x) == len(pst.eti_pdf) - assert len(pst.hdi_x) == len(pst.hdi_pdf) - assert tuple(pst.eti_x[[0, -1]]) == (0, 20) - assert tuple(pst.hdi_x[[0, -1]]) == (0, 20) - assert (numpy.isclose(scipy.integrate.cumtrapz(pst.eti_pdf, pst.eti_x)[-1], 1, atol=0.0001)) - assert (numpy.isclose(scipy.integrate.cumtrapz(pst.hdi_pdf, pst.hdi_x)[-1], 1, atol=0.0001)) - assert pst.eti_lower == pst.hdi_lower == 0 - assert pst.eti_upper == pst.hdi_upper == 20 - assert pst.eti_prob == 1 - assert pst.hdi_prob == 1 - - # check trimming to [2.5,97.5] interval - pst = em.infer_independent(y=[1, 2], lower=0, upper=20, steps=1775, ci_prob=0.95) - - assert len(pst.eti_x) == len(pst.eti_pdf) - assert len(pst.hdi_x) == len(pst.hdi_pdf) - assert numpy.isclose(pst.eti_prob, 0.95, atol=0.0001) - assert numpy.isclose(pst.hdi_prob, 0.95, atol=0.0001) - assert numpy.isclose(scipy.integrate.cumtrapz(pst.eti_pdf, pst.eti_x)[-1], 0.95, atol=0.0001) - assert numpy.isclose(scipy.integrate.cumtrapz(pst.hdi_pdf, pst.hdi_x)[-1], 0.95, atol=0.0001) - assert pst.eti_lower == pst.eti_x[0] - assert pst.eti_upper == pst.eti_x[-1] - assert pst.hdi_lower == pst.hdi_x[0] - assert pst.hdi_upper == pst.hdi_x[-1] - - # check that error are raised by wrong input - with pytest.raises(ValueError): - _ = em.infer_independent(y=1, lower=0, upper=20, steps=1000, ci_prob=(-1)) - with pytest.raises(ValueError): - _ = em.infer_independent(y=1, lower=0, upper=20, steps=1000, ci_prob=(97.5)) - pass - - @pytest.mark.skipif(not HAS_PYMC, reason='requires PyMC') - def test_symbolic_loglikelihood_checks_and_warnings(self): - cmodel = _TestPolynomialModel(independent_key='S', dependent_key='A', mu_degree=1, scale_degree=1) - cmodel.theta_fitted = [0, 1, 0.1, 1, 5] - - # create test data - x_true = numpy.array([1,2,3,4,5]) - y_obs = cmodel.predict_dependent(x_true)[0] - - with pm.Model(): - x_hat = pm.Uniform("x_hat", shape=5) - with pytest.raises(ValueError, match="`name` must be specified"): - cmodel.loglikelihood(x=x_hat, y=y_obs) - - with pytest.warns(DeprecationWarning, match="Use `name` instead"): - cmodel.loglikelihood(x=x_hat, y=y_obs, replicate_id='A01', dependent_key='A') - pass - - @pytest.mark.skipif(not HAS_PYMC, reason='requires PyMC') - def test_symbolic_loglikelihood(self): - cmodel = _TestPolynomialModel(independent_key='S', dependent_key='A', mu_degree=1, scale_degree=1) - cmodel.theta_fitted = [0, 1, 0.1, 1, 3] - - # create test data - x_true = numpy.array([1,2,3,4,5]) - y_obs = cmodel.predict_dependent(x_true)[0] - - x_hat = at.vector() - x_hat.tag.test_value = x_true - L = cmodel.loglikelihood(x=x_hat, y=y_obs, name='L_A01_A') - assert isinstance(L, at.TensorVariable) - assert L.ndim == 0 - - # compare the two loglikelihood computation methods - x_test = numpy.random.normal(x_true, scale=0.1) - actual = L.eval({ - x_hat: x_test - }) - expected = cmodel.loglikelihood(x=x_test, y=y_obs) - assert numpy.ndim(expected) == 0 - assert numpy.ndim(actual) == 0 - numpy.testing.assert_almost_equal(actual, expected, 6) - pass - - @pytest.mark.skipif(not HAS_PYMC, reason='requires PyMC') - def test_symbolic_loglikelihood_in_modelcontext(self): - cmodel = _TestPolynomialModel(independent_key='S', dependent_key='A', mu_degree=1, scale_degree=1) - cmodel.theta_fitted = [0, 0.5, 0.1, 1, 7] - - # create test data - x_true = numpy.array([1,2,3,4,5]) - y_obs = cmodel.predict_dependent(x_true)[0] - - # create a PyMC model using the calibration model - with pm.Model() as pmodel: - x_hat = pm.Uniform("x_hat", 0, 1, shape=x_true.shape, transform=None) - L = cmodel.loglikelihood(x=x_hat, y=y_obs, name='L_A01_A') - assert isinstance(L, at.TensorVariable) - assert L.ndim == 0 - - # PyMC v4 returns the RV, but for .eval() we need the RV-value-variable - if pm.__version__[0] != "3": - x_hat = pmodel.rvs_to_values[x_hat] - - # compare the two loglikelihood computation methods - x_test = numpy.random.normal(x_true, scale=0.1) - actual = L.eval({ - x_hat: x_test - }) - expected = cmodel.loglikelihood(x=x_test, y=y_obs) - assert numpy.ndim(expected) == 0 - assert numpy.ndim(actual) == 0 - numpy.testing.assert_almost_equal(actual, expected, 6) - pass - - @pytest.mark.parametrize("x", [ - numpy.array([1,2,3]), - 4, - ]) - @pytest.mark.parametrize("y", [ - numpy.array([2,4,8]), - 5, - ]) - def test_loglikelihood(self, x, y): - cmodel = _TestPolynomialModel(independent_key='S', dependent_key='OD', mu_degree=1, scale_degree=1) - cmodel.theta_fitted = [0, 1, 0.1, 1.6, 3] - - actual = cmodel.loglikelihood(y=y, x=x) - assert numpy.ndim(actual) == 0 - mu, scale, df = cmodel.predict_dependent(x, theta=cmodel.theta_fitted) - expected = numpy.sum(stats.t.logpdf(x=y, loc=mu, scale=scale, df=df)) - numpy.testing.assert_equal(actual, expected) - return - - def test_loglikelihood_exceptions(self): - cmodel = _TestPolynomialModel(independent_key='S', dependent_key='OD', mu_degree=1, scale_degree=1) - with pytest.raises(Exception, match="No parameter vector"): - cmodel.loglikelihood(y=[2,3], x=[4,5]) - - cmodel.theta_fitted = [0, 1, 0.1, 1.6, 2] - - with pytest.raises(TypeError): - cmodel.loglikelihood(4, x=[2,3]) - with pytest.raises(ValueError, match="Input x must be"): - cmodel.loglikelihood(y=[2,3], x="hello") - with pytest.raises(ValueError, match="operands could not be broadcast"): - cmodel.loglikelihood(y=[1,2,3], x=[1,2]) - return - - def test_likelihood(self): - # use a linear model with intercept 1 and slope 0.5 - cmodel = _TestPolynomialModel(independent_key="I", dependent_key="D", mu_degree=1) - cmodel.theta_fitted = [1, 0.5, 0.5, 4] - - assert numpy.isscalar(cmodel.likelihood(y=2, x=3)) - assert numpy.isscalar(cmodel.likelihood(y=[2, 3], x=[3, 4])) - - with pytest.raises(ValueError, match="operands could not be broadcast"): - cmodel.likelihood(y=[1,2,3], x=[1,2]) - - x_dense = numpy.linspace(0, 4, 501) - actual = cmodel.likelihood(x=x_dense, y=2, scan_x=True) - assert numpy.ndim(actual) == 1 - # the maximum likelihood should be at x=2 - assert x_dense[numpy.argmax(actual)] == 2 - pass - - @pytest.mark.xfail(reason="Draft test case, see https://github.com/JuBiotech/calibr8/issues/15") - def test_likelihood_nobroadcasting_fallback(self): - class _TestSwitchableBroadcastingModel(calibr8.CalibrationModel, calibr8.NormalNoise): - def __init__(self): - self.x_shapes = [] - self.can_broadcast = False - super().__init__(independent_key="I", dependent_key="D", theta_names="i,s,sd".split(","), ndim=1) - - def predict_dependent(self, x, *, theta=None): - if theta is None: - theta = self.theta_fitted - i, s, sd = self.theta_fitted - # This part broadcasts just fine - x = numpy.array(x) - mu = i + x * s - return mu, sd - - def loglikelihood(self, *, y, x, **kwargs): - # This overrides the native CalibrationModel.loglikelihood with one - # that can be externally set to return non-broadcasted results. - LL_broadcasted = super().loglikelihood(y=y, x=x, **kwargs) - if not self.can_broadcast: - return LL_broadcasted.sum() - return LL_broadcasted - - cm = _TestSwitchableBroadcastingModel() - cm.theta_fitted = [0.5, 0.6, 0.7] - - X = numpy.array([1, 2, 3]) - Y = [1.1, 1.7, 2.3] - y_obs = [1.0, 1.7, 2.4] - - # Check the prediction of the test model - for x, y in zip(X, Y): - assert cm.predict_dependent(x) == (y, 0.7) - - # The predict_dependent can broadcast - mu, sd = cm.predict_dependent(X) - numpy.testing.assert_array_equal(mu, Y) - assert sd == 0.7 - - # Test the switching between can_broadcast modes - cm.can_broadcast = True - LL = cm.loglikelihood(x=X[..., None], y=y_obs) - assert numpy.shape(LL) == (3,) - - cm.can_broadcast = False - LL = cm.loglikelihood(x=X[..., None], y=y_obs) - assert numpy.shape(LL) == () - - # The CalibrationModel.likelihood should give the same results either way. - cm.can_broadcast = True - L_broadcasted = cm.likelihood(x=X[..., None], y=y_obs, scan_x=True) - cm.can_broadcast = False - L_looped = cm.likelihood(x=X[..., None], y=y_obs, scan_x=True) - - numpy.testing.assert_array_equal(L_broadcasted, L_looped) - - # Of course the values should also be correct. - L_expected = numpy.exp([ - cm.scipy_dist.logpdf(x=y_obs, loc=mui, scale=sd).sum() - for mui in mu - ]) - assert numpy.shape(L_expected) == (3,) - numpy.testing.assert_array_equal(L_broadcasted, L_expected) - pass - - -class TestBaseAsymmetricLogisticModelT: - def test_predict_dependent(self): - x = numpy.array([1,2,3]) - theta = [0, 4, 2, 1, 1, 0, 2, 1.4] - cmodel = _TestLogisticModel(independent_key='S', dependent_key='OD', scale_degree=1) - cmodel.theta_fitted = theta - with pytest.raises(TypeError): - _ = cmodel.predict_dependent(x, theta) - mu, scale, df = cmodel.predict_dependent(x) - numpy.testing.assert_array_equal(mu, calibr8.asymmetric_logistic(x, theta)) - numpy.testing.assert_array_equal(scale, 0 + 2 * mu) - assert df == 1.4 - return - - def test_predict_independent(self): - cmodel = _TestLogisticModel(independent_key='S', dependent_key='OD', scale_degree=1) - cmodel.theta_fitted = [0, 4, 2, 1, 1, 2, 1.43, 5] - x_original = numpy.array([4, 5, 6]) - mu, sd, df = cmodel.predict_dependent(x_original) - x_predicted = cmodel.predict_independent(y=mu) - assert (numpy.allclose(x_predicted, x_original)) - return - - class TestOptimization: def _get_test_model(self): theta_mu = (0.5, 1.4) @@ -1330,3 +1225,146 @@ def test_fit_pygmo(self, caplog): numpy.testing.assert_array_equal(em.cal_independent, x) numpy.testing.assert_array_equal(em.cal_dependent, y) pass + + +class TestContribBase: + def test_cant_instantiate_base_models(self): + with pytest.raises(TypeError, match=r"constructor must not have any required \(kw\)arguments"): + calibr8.BaseModelT(independent_key='I', dependent_key='D', theta_names=["a", "b"]) + with pytest.raises(TypeError, match=r"constructor must not have any required \(kw\)arguments"): + calibr8.BaseAsymmetricLogisticT(independent_key='I', dependent_key='D') + with pytest.raises(TypeError, match=r"constructor must not have any required \(kw\)arguments"): + calibr8.BasePolynomialModelT(independent_key='I', dependent_key='D', mu_degree=1, scale_degree=1) + pass + + +class TestBasePolynomialModelT: + def test_exceptions(self): + with pytest.raises(ValueError, match="are useless"): + _TestPolynomialModel(independent_key='I', dependent_key='D', mu_degree=0) + pass + + @pytest.mark.parametrize("mu_degree,scale_degree", [(1, 0), (1, 1), (1, 2), (2, 0)]) + def test_predict_dependent(self, mu_degree, scale_degree): + theta_mu = (2.2, 1.2, 0.2)[:mu_degree+1] + theta_scale = (3.1, 0.4, 0.2)[:scale_degree+1] + theta = theta_mu + theta_scale + (1,) + + em = _TestPolynomialModel(independent_key='I', dependent_key='D', mu_degree=mu_degree, scale_degree=scale_degree) + assert len(em.theta_names) == mu_degree+1 + scale_degree+1 + 1 + assert len(em.theta_names) == len(theta) + + x = numpy.linspace(0, 10, 3) + mu, scale, df = em.predict_dependent(x, theta=theta) + + + expected = numpy.polyval(theta_mu[::-1], x) + numpy.testing.assert_array_equal(mu, expected) + + expected = numpy.polyval(theta_scale[::-1], mu) + numpy.testing.assert_array_equal(scale, expected) + + numpy.testing.assert_array_equal(df, 1) + pass + + @pytest.mark.parametrize("mu_degree,scale_degree", [(1, 0), (1, 1), (1, 2), (2, 0)]) + def test_predict_independent(self, mu_degree, scale_degree): + theta_mu = (2.2, 1.2, 0.2)[:mu_degree+1] + theta_scale = (3.1, 0.4, 0.2)[:scale_degree+1] + theta = theta_mu + theta_scale + (1,) + + em = _TestPolynomialModel(independent_key='I', dependent_key='D', mu_degree=mu_degree, scale_degree=scale_degree) + em.theta_fitted = theta + assert len(em.theta_names) == mu_degree+1 + scale_degree+1 + 1 + assert len(em.theta_names) == len(theta) + + x = numpy.linspace(0, 10, 7) + mu, _, _ = em.predict_dependent(x, theta=theta) + + if mu_degree < 2: + x_inverse = em.predict_independent(mu) + numpy.testing.assert_array_almost_equal(x_inverse, x) + else: + with pytest.raises(NotImplementedError, match="higher order polynomials"): + em.predict_independent(mu) + pass + + +class TestBaseAsymmetricLogisticModelT: + @pytest.mark.parametrize("scale_degree", [0, 1, 2]) + def test_predict_dependent(self, scale_degree): + theta_mu = (-0.5, 0.5, 1, 1, -1) + theta_scale = (3.1, 0.4, 0.2)[:scale_degree+1] + theta = theta_mu + theta_scale + (1,) + + em = _TestLogisticModel(independent_key='I', dependent_key='D', scale_degree=scale_degree) + assert len(em.theta_names) == 5 + scale_degree+1 + 1 + assert len(em.theta_names) == len(theta) + + x = numpy.linspace(1, 5, 3) + mu, scale, df = em.predict_dependent(x, theta=theta) + + expected = calibr8.asymmetric_logistic(x, theta_mu) + numpy.testing.assert_array_equal(mu, expected) + + expected = numpy.polyval(theta_scale[::-1], mu) + numpy.testing.assert_array_equal(scale, expected) + + numpy.testing.assert_array_equal(df, 1) + pass + + @pytest.mark.parametrize("scale_degree", [0, 1, 2]) + def test_predict_independent(self, scale_degree): + theta_mu = (-0.5, 0.5, 1, 1, -1) + theta_scale = (3.1, 0.4, 0.2)[:scale_degree+1] + theta = theta_mu + theta_scale + (1,) + + em = _TestLogisticModel(independent_key='I', dependent_key='D', scale_degree=scale_degree) + em.theta_fitted = theta + + x = numpy.linspace(1, 5, 7) + mu, _, _ = em.predict_dependent(x) + + x_inverse = em.predict_independent(mu) + numpy.testing.assert_array_almost_equal(x_inverse, x) + pass + + +class TestBaseLogIndependentAsymmetricModelT: + @pytest.mark.parametrize("scale_degree", [0, 1, 2]) + def test_predict_dependent(self, scale_degree): + theta_mu = (-0.5, 0.5, 1, 1, -1) + theta_scale = (3.1, 0.4, 0.2)[:scale_degree+1] + theta = theta_mu + theta_scale + (1,) + + em = _TestLogIndependentLogisticModel(independent_key='I', dependent_key='D', scale_degree=scale_degree) + assert len(em.theta_names) == 5 + scale_degree+1 + 1 + assert len(em.theta_names) == len(theta) + + x = numpy.linspace(1, 5, 3) + mu, scale, df = em.predict_dependent(x, theta=theta) + + expected = calibr8.xlog_asymmetric_logistic(x, theta_mu) + numpy.testing.assert_array_equal(mu, expected) + + expected = numpy.polyval(theta_scale[::-1], mu) + numpy.testing.assert_array_equal(scale, expected) + + numpy.testing.assert_array_equal(df, 1) + pass + + @pytest.mark.parametrize("scale_degree", [0, 1, 2]) + def test_predict_independent(self, scale_degree): + theta_mu = (-0.5, 0.5, 1, 1, -1) + theta_scale = (3.1, 0.4, 0.2)[:scale_degree+1] + theta = theta_mu + theta_scale + (1,) + + em = _TestLogIndependentLogisticModel(independent_key='I', dependent_key='D', scale_degree=scale_degree) + em.theta_fitted = theta + + x = numpy.linspace(1, 5, 7) + mu, _, _ = em.predict_dependent(x) + + x_inverse = em.predict_independent(mu) + numpy.testing.assert_array_almost_equal(x_inverse, x) + pass