From 2848fe4e608e43492b03ae1ef2bf6a0479be1044 Mon Sep 17 00:00:00 2001 From: Erik Bernhardsson Date: Mon, 22 Jul 2019 22:53:13 -0400 Subject: [PATCH 1/6] move all params to sit on the constructor of each class --- convoys/regression.py | 43 ++++++++++++++++++++++++------------------- 1 file changed, 24 insertions(+), 19 deletions(-) diff --git a/convoys/regression.py b/convoys/regression.py index 7198af4..095851c 100644 --- a/convoys/regression.py +++ b/convoys/regression.py @@ -147,10 +147,12 @@ class GeneralizedGamma(RegressionModel): to sample from the full posterior in order to generate uncertainty estimates for all parameters. ''' - def __init__(self, ci=False): + def __init__(self, ci=False, fix_k=None, fix_p=None): self._ci = ci + self._fix_k = fix_k + self._fix_p = fix_p - def fit(self, X, B, T, W=None, fix_k=None, fix_p=None): + def fit(self, X, B, T, W=None): '''Fits the model. :param X: numpy matrix of shape :math:`k \\cdot n` @@ -176,9 +178,9 @@ def fit(self, X, B, T, W=None, fix_k=None, fix_p=None): # a, b, alpha_1...alpha_k, beta_1...beta_k) # Generalized Gamma is a bit sensitive to the starting point! x0 = numpy.zeros(6+2*n_features) - x0[0] = +1 if fix_k is None else log(fix_k) - x0[1] = -1 if fix_p is None else log(fix_p) - args = (X, B, T, W, fix_k, fix_p, True) + x0[0] = +1 if self._fix_k is None else log(self._fix_k) + x0[1] = -1 if self._fix_p is None else log(self._fix_p) + args = (X, B, T, W, self._fix_k, self._fix_p, True) # Callback for progress to stdout sys.stdout.write('\n') @@ -203,10 +205,10 @@ def callback(LL, value_history=[]): result = {'map': res.x} # TODO: should not use fixed k/p as search parameters - if fix_k: - result['map'][0] = log(fix_k) - if fix_p: - result['map'][1] = log(fix_p) + if self._fix_k: + result['map'][0] = log(self._fix_k) + if self._fix_p: + result['map'][1] = log(self._fix_p) # Make sure we're in a local minimum gradient = jac(result['map']) @@ -238,10 +240,10 @@ def callback(LL, value_history=[]): sys.stdout.write('\n') result['samples'] = sampler.chain[:, n_burnin:, :] \ .reshape((-1, dim)).T - if fix_k: - result['samples'][0, :] = log(fix_k) - if fix_p: - result['samples'][1, :] = log(fix_p) + if self._fix_k: + result['samples'][0, :] = log(self._fix_k) + if self._fix_p: + result['samples'][1, :] = log(self._fix_p) self.params = {k: { 'k': exp(data[0]), @@ -324,8 +326,9 @@ class Exponential(GeneralizedGamma): from the initial state to converted or dead is constant. See documentation for :class:`GeneralizedGamma`.''' - def fit(self, X, B, T, W=None): - super(Exponential, self).fit(X, B, T, W, fix_k=1, fix_p=1) + def __init__(self, *args, **kwargs): + kwargs.update(dict(fix_k=1, fix_p=1)) + super(Exponential, self).__init__(*args, **kwargs) class Weibull(GeneralizedGamma): @@ -340,8 +343,9 @@ class Weibull(GeneralizedGamma): :math:`f(t) = p\\lambda(t\\lambda)^{p-1}\\exp(-(t\\lambda)^p)` See documentation for :class:`GeneralizedGamma`.''' - def fit(self, X, B, T, W=None): - super(Weibull, self).fit(X, B, T, W, fix_k=1) + def __init__(self, *args, **kwargs): + kwargs.update(dict(fix_k=1)) + super(Weibull, self).__init__(*args, **kwargs) class Gamma(GeneralizedGamma): @@ -359,5 +363,6 @@ class Gamma(GeneralizedGamma): :math:`f(t) = \\lambda^k t^{k-1} \\exp(-x\\lambda) / \\Gamma(k)` See documentation for :class:`GeneralizedGamma`.''' - def fit(self, X, B, T, W=None): - super(Gamma, self).fit(X, B, T, W, fix_p=1) + def __init__(self, *args, **kwargs): + kwargs.update(dict(fix_p=1)) + super(Gamma, self).__init__(*args, **kwargs) From ff4581254b3c4bdc7b695592e4f022f108010088 Mon Sep 17 00:00:00 2001 From: Erik Bernhardsson Date: Mon, 22 Jul 2019 22:55:25 -0400 Subject: [PATCH 2/6] make hierarchical an option --- convoys/regression.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/convoys/regression.py b/convoys/regression.py index 095851c..2bd6c49 100644 --- a/convoys/regression.py +++ b/convoys/regression.py @@ -147,10 +147,11 @@ class GeneralizedGamma(RegressionModel): to sample from the full posterior in order to generate uncertainty estimates for all parameters. ''' - def __init__(self, ci=False, fix_k=None, fix_p=None): + def __init__(self, ci=False, fix_k=None, fix_p=None, hierarchical=True): self._ci = ci self._fix_k = fix_k self._fix_p = fix_p + self._hierarchical = hierarchical def fit(self, X, B, T, W=None): '''Fits the model. @@ -180,7 +181,7 @@ def fit(self, X, B, T, W=None): x0 = numpy.zeros(6+2*n_features) x0[0] = +1 if self._fix_k is None else log(self._fix_k) x0[1] = -1 if self._fix_p is None else log(self._fix_p) - args = (X, B, T, W, self._fix_k, self._fix_p, True) + args = (X, B, T, W, self._fix_k, self._fix_p, self._hierarchical) # Callback for progress to stdout sys.stdout.write('\n') From 386baf138300292f31e5e92175ac067c9769457f Mon Sep 17 00:00:00 2001 From: Erik Bernhardsson Date: Mon, 22 Jul 2019 23:28:56 -0400 Subject: [PATCH 3/6] add linear model --- convoys/regression.py | 28 ++++++++++++++++++---------- test_convoys.py | 6 ++++++ 2 files changed, 24 insertions(+), 10 deletions(-) diff --git a/convoys/regression.py b/convoys/regression.py index 2bd6c49..7dc3f4e 100644 --- a/convoys/regression.py +++ b/convoys/regression.py @@ -17,8 +17,8 @@ 'GeneralizedGamma'] -def generalized_gamma_LL(x, X, B, T, W, fix_k, fix_p, - hierarchical, callback=None): +def generalized_gamma_loss(x, X, B, T, W, fix_k, fix_p, + hierarchical, flavor, callback=None): k = exp(x[0]) if fix_k is None else fix_k p = exp(x[1]) if fix_p is None else fix_p log_sigma_alpha = x[2] @@ -29,15 +29,20 @@ def generalized_gamma_LL(x, X, B, T, W, fix_k, fix_p, alpha = x[6:6+n_features] beta = x[6+n_features:6+2*n_features] lambd = exp(dot(X, alpha)+a) - c = expit(dot(X, beta)+b) # PDF: p*lambda^(k*p) / gamma(k) * t^(k*p-1) * exp(-(x*lambda)^p) log_pdf = log(p) + (k*p) * log(lambd) - gammaln(k) \ + (k*p-1) * log(T) - (T*lambd)**p cdf = gammainc(k, (T*lambd)**p) - LL_observed = log(c) + log_pdf - LL_censored = log((1-c) + c * (1 - cdf)) + if flavor == 'logistic': # Log-likelihood with sigmoid + c = expit(dot(X, beta)+b) + LL_observed = log(c) + log_pdf + LL_censored = log((1 - c) + c * (1 - cdf)) + elif flavor == 'linear': # L2 loss, linear + c = dot(X, beta)+b + LL_observed = -(1 - c)**2 + log_pdf + LL_censored = -c**2 * cdf - (1 - c)**2 * (1 - cdf) LL_data = sum( W * B * LL_observed + @@ -147,11 +152,13 @@ class GeneralizedGamma(RegressionModel): to sample from the full posterior in order to generate uncertainty estimates for all parameters. ''' - def __init__(self, ci=False, fix_k=None, fix_p=None, hierarchical=True): + def __init__(self, ci=False, fix_k=None, fix_p=None, hierarchical=True, + flavor='logistic'): self._ci = ci self._fix_k = fix_k self._fix_p = fix_p self._hierarchical = hierarchical + self._flavor = flavor def fit(self, X, B, T, W=None): '''Fits the model. @@ -181,7 +188,8 @@ def fit(self, X, B, T, W=None): x0 = numpy.zeros(6+2*n_features) x0[0] = +1 if self._fix_k is None else log(self._fix_k) x0[1] = -1 if self._fix_p is None else log(self._fix_p) - args = (X, B, T, W, self._fix_k, self._fix_p, self._hierarchical) + args = (X, B, T, W, self._fix_k, self._fix_p, + self._hierarchical, self._flavor) # Callback for progress to stdout sys.stdout.write('\n') @@ -193,8 +201,8 @@ def callback(LL, value_history=[]): sys.stdout.flush() # Define objective and use automatic differentiation - f = lambda x: -generalized_gamma_LL(x, *args, callback=callback) - jac = autograd.grad(lambda x: -generalized_gamma_LL(x, *args)) + f = lambda x: -generalized_gamma_loss(x, *args, callback=callback) + jac = autograd.grad(lambda x: -generalized_gamma_loss(x, *args)) # Find the maximum a posteriori of the distribution res = scipy.optimize.minimize(f, x0, jac=jac, method='SLSQP', @@ -225,7 +233,7 @@ def callback(LL, value_history=[]): sampler = emcee.EnsembleSampler( nwalkers=n_walkers, dim=dim, - lnpostfn=generalized_gamma_LL, + lnpostfn=generalized_gamma_loss, args=args, ) mcmc_initial_noise = 1e-3 diff --git a/test_convoys.py b/test_convoys.py index b140033..bc39824 100644 --- a/test_convoys.py +++ b/test_convoys.py @@ -94,6 +94,12 @@ def test_exponential_regression_model(c=0.3, lambd=0.1, n=10000): assert model.cdf([1], 0).shape == () assert model.cdf([1], [0, 1, 2, 3]).shape == (4,) + # Fit a linear model + model = convoys.regression.Exponential(ci=False, flavor='linear') + model.fit(X, B, T) + model_c = model.params['map']['b'] + model.params['map']['beta'][0] + assert 0.9*c < model_c < 1.1*c + @flaky.flaky def test_weibull_regression_model(cs=[0.3, 0.5, 0.7], From d9c229f0e43df2ea92d5df443d6619e14fe03b93 Mon Sep 17 00:00:00 2001 From: Erik Bernhardsson Date: Mon, 22 Jul 2019 23:49:49 -0400 Subject: [PATCH 4/6] think ll makes more sense now, but added a breaking weibull test --- convoys/regression.py | 2 +- test_convoys.py | 7 +++++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/convoys/regression.py b/convoys/regression.py index 7dc3f4e..e82b008 100644 --- a/convoys/regression.py +++ b/convoys/regression.py @@ -42,7 +42,7 @@ def generalized_gamma_loss(x, X, B, T, W, fix_k, fix_p, elif flavor == 'linear': # L2 loss, linear c = dot(X, beta)+b LL_observed = -(1 - c)**2 + log_pdf - LL_censored = -c**2 * cdf - (1 - c)**2 * (1 - cdf) + LL_censored = -c**2 * cdf LL_data = sum( W * B * LL_observed + diff --git a/test_convoys.py b/test_convoys.py index bc39824..dc432ac 100644 --- a/test_convoys.py +++ b/test_convoys.py @@ -127,6 +127,13 @@ def test_weibull_regression_model(cs=[0.3, 0.5, 0.7], x = [int(r == j) for j in range(len(cs))] assert 0.80 * c < model.cdf(x, float('inf')) < 1.30 * c + # Fit a linear model + model = convoys.regression.Weibull(ci=False, flavor='linear') + model.fit(X, B, T) + model_cs = model.params['map']['b'] + model.params['map']['beta'] + for model_c, c in zip(model_cs, cs): + assert 0.8 * c < model_c < 1.2 * c + @flaky.flaky def test_gamma_regression_model(c=0.3, lambd=0.1, k=3.0, n=10000): From 5e34d8e0e8cf4d6a09265cd791c94dd67fe80de3 Mon Sep 17 00:00:00 2001 From: Erik Bernhardsson Date: Thu, 1 Aug 2019 17:27:11 -0400 Subject: [PATCH 5/6] add a more complicated test for linear regression --- convoys/regression.py | 2 +- test_convoys.py | 27 +++++++++++++++++++++++++++ 2 files changed, 28 insertions(+), 1 deletion(-) diff --git a/convoys/regression.py b/convoys/regression.py index e82b008..6958ebb 100644 --- a/convoys/regression.py +++ b/convoys/regression.py @@ -42,7 +42,7 @@ def generalized_gamma_loss(x, X, B, T, W, fix_k, fix_p, elif flavor == 'linear': # L2 loss, linear c = dot(X, beta)+b LL_observed = -(1 - c)**2 + log_pdf - LL_censored = -c**2 * cdf + LL_censored = -(c*cdf)**2 LL_data = sum( W * B * LL_observed + diff --git a/test_convoys.py b/test_convoys.py index dc432ac..2bce37a 100644 --- a/test_convoys.py +++ b/test_convoys.py @@ -149,6 +149,33 @@ def test_gamma_regression_model(c=0.3, lambd=0.1, k=3.0, n=10000): assert 0.80*c < model.cdf([1], float('inf')) < 1.30*c assert 0.80*k < numpy.mean(model.params['map']['k']) < 1.30*k + # Fit a linear model + model = convoys.regression.Gamma(ci=False, flavor='linear') + model.fit(X, B, T) + model_c = model.params['map']['b'] + model.params['map']['beta'][0] + assert 0.9*c < model_c < 1.1*c + + +@flaky.flaky +def test_linear_model(n=10000, k=5, lambd=0.1): + # Generate data with little censoring + # The coefficients should be quite close to their actual value + cs = numpy.random.dirichlet(numpy.ones(k)) + X = numpy.random.binomial(n=1, p=0.5, size=(n, k)) + C = numpy.random.rand(n) < numpy.dot(X, cs.T) + N = scipy.stats.uniform.rvs(scale=20./lambd, size=(n,)) + E = numpy.array([sample_weibull(k, lambd) for r in range(n)]) + B, T = generate_censored_data(N, E, C) + + model = convoys.regression.Weibull(ci=False, flavor='linear') + model.fit(X, B, T) + model_cs = model.params['map']['b'] + model.params['map']['beta'] + for model_c, c in zip(model_cs, cs): + assert c - 0.03 < model_c < c + 0.03 + model_lambds = numpy.exp(model.params['map']['a'] + model.params['map']['alpha']) + for model_lambd in model_lambds: + assert 0.97*lambd < model_lambd < 1.03*lambd + @flaky.flaky def test_exponential_pooling(c=0.5, lambd=0.01, n=10000, ks=[1, 2, 3]): From 2c9334b20721705e57fa8d215d9e4f2982f96fd6 Mon Sep 17 00:00:00 2001 From: Erik Bernhardsson Date: Thu, 1 Aug 2019 17:29:01 -0400 Subject: [PATCH 6/6] hound --- test_convoys.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test_convoys.py b/test_convoys.py index 2bce37a..63eb491 100644 --- a/test_convoys.py +++ b/test_convoys.py @@ -172,7 +172,8 @@ def test_linear_model(n=10000, k=5, lambd=0.1): model_cs = model.params['map']['b'] + model.params['map']['beta'] for model_c, c in zip(model_cs, cs): assert c - 0.03 < model_c < c + 0.03 - model_lambds = numpy.exp(model.params['map']['a'] + model.params['map']['alpha']) + model_lambds = numpy.exp(model.params['map']['a'] + + model.params['map']['alpha']) for model_lambd in model_lambds: assert 0.97*lambd < model_lambd < 1.03*lambd