diff --git a/GPy/core/domains.py b/GPy/core/domains.py new file mode 100644 index 00000000..dfc880f6 --- /dev/null +++ b/GPy/core/domains.py @@ -0,0 +1,10 @@ +''' +Created on 4 Jun 2013 + +@author: maxz +''' + +REAL = 'real' +POSITIVE = "positive" +NEGATIVE = 'negative' +BOUNDED = 'bounded' diff --git a/GPy/core/model.py b/GPy/core/model.py index 54b1fdd6..3a1d08e9 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -41,16 +41,16 @@ class model(parameterised): Arguments --------- regexp -- string, regexp, or integer array - what -- instance of a prior class + what -- instance of a Prior class Notes ----- - Asserts that the prior is suitable for the constraint. If the + Asserts that the Prior is suitable for the constraint. If the wrong constraint is in place, an error is raised. If no constraint is in place, one is added (warning printed). - For tied parameters, the prior will only be "counted" once, thus - a prior object is only inserted on the first tied index + For tied parameters, the Prior will only be "counted" once, thus + a Prior object is only inserted on the first tied index """ if self.priors is None: self.priors = [None for i in range(self._get_params().size)] @@ -60,24 +60,24 @@ class model(parameterised): # check tied situation tie_partial_matches = [tie for tie in self.tied_indices if (not set(tie).isdisjoint(set(which))) & (not set(tie) == set(which))] if len(tie_partial_matches): - raise ValueError, "cannot place prior across partial ties" + raise ValueError, "cannot place Prior across partial ties" tie_matches = [tie for tie in self.tied_indices if set(which) == set(tie) ] if len(tie_matches) > 1: - raise ValueError, "cannot place prior across multiple ties" + raise ValueError, "cannot place Prior across multiple ties" elif len(tie_matches) == 1: - which = which[:1] # just place a prior object on the first parameter + which = which[:1] # just place a Prior object on the first parameter # check constraints are okay if what.domain is POSITIVE: - constrained_positive_indices = [i for i, t in zip(self.constrained_indices, self.constraints) if t.domain == POSITIVE] + constrained_positive_indices = [i for i, t in zip(self.constrained_indices, self.constraints) if t.domain is POSITIVE] if len(constrained_positive_indices): constrained_positive_indices = np.hstack(constrained_positive_indices) else: constrained_positive_indices = np.zeros(shape=(0,)) bad_constraints = np.setdiff1d(self.all_constrained_indices(), constrained_positive_indices) - assert not np.any(which[:, None] == bad_constraints), "constraint and prior incompatible" + assert not np.any(which[:, None] == bad_constraints), "constraint and Prior incompatible" unconst = np.setdiff1d(which, constrained_positive_indices) if len(unconst): print "Warning: constraining parameters to be positive:" @@ -85,11 +85,11 @@ class model(parameterised): print '\n' self.constrain_positive(unconst) elif what.domain is REAL: - assert not np.any(which[:, None] == self.all_constrained_indices()), "constraint and prior incompatible" + assert not np.any(which[:, None] == self.all_constrained_indices()), "constraint and Prior incompatible" else: - raise ValueError, "prior not recognised" + raise ValueError, "Prior not recognised" - # store the prior in a local list + # store the Prior in a local list for w in which: self.priors[w] = what @@ -136,13 +136,13 @@ class model(parameterised): def randomize(self): """ Randomize the model. - Make this draw from the prior if one exists, else draw from N(0,1) + Make this draw from the Prior if one exists, else draw from N(0,1) """ # first take care of all parameters (from N(0,1)) x = self._get_params_transformed() x = np.random.randn(x.size) self._set_params_transformed(x) - # now draw from prior where possible + # now draw from Prior where possible x = self._get_params() if self.priors is not None: [np.put(x, i, p.rvs(1)) for i, p in enumerate(self.priors) if not p is None] @@ -284,7 +284,7 @@ class model(parameterised): def Laplace_covariance(self): """return the covariance matric of a Laplace approximatino at the current (stationary) point""" - # TODO add in the prior contributions for MAP estimation + # TODO add in the Prior contributions for MAP estimation # TODO fix the hessian for tied, constrained and fixed components if hasattr(self, 'log_likelihood_hessian'): A = -self.log_likelihood_hessian() @@ -323,14 +323,14 @@ class model(parameterised): log_prior = self.log_prior() obj_funct = '\nLog-likelihood: {0:.3e}'.format(log_like) if len(''.join(strs)) != 0: - obj_funct += ', Log prior: {0:.3e}, LL+prior = {0:.3e}'.format(log_prior, log_like + log_prior) + obj_funct += ', Log Prior: {0:.3e}, LL+Prior = {0:.3e}'.format(log_prior, log_like + log_prior) obj_funct += '\n\n' s[0] = obj_funct + s[0] s[0] += "|{h:^{col}}".format(h='Prior', col=width) s[1] += '-' * (width + 1) for p in range(2, len(strs) + 2): - s[p] += '|{prior:^{width}}'.format(prior=strs[p - 2], width=width) + s[p] += '|{Prior:^{width}}'.format(Prior=strs[p - 2], width=width) return '\n'.join(s) diff --git a/GPy/core/priors.py b/GPy/core/priors.py index 74ca63bf..b86024d0 100644 --- a/GPy/core/priors.py +++ b/GPy/core/priors.py @@ -9,7 +9,7 @@ from ..util.linalg import pdinv from GPy.core.domains import REAL, POSITIVE import warnings -class prior: +class Prior: domain = None def pdf(self, x): return np.exp(self.lnpdf(x)) @@ -22,7 +22,7 @@ class prior: pb.plot(xx, self.pdf(xx), 'r', linewidth=2) -class Gaussian(prior): +class Gaussian(Prior): """ Implementation of the univariate Gaussian probability function, coupled with random variables. @@ -52,7 +52,7 @@ class Gaussian(prior): return np.random.randn(n) * self.sigma + self.mu -class log_Gaussian(prior): +class LogGaussian(Prior): """ Implementation of the univariate *log*-Gaussian probability function, coupled with random variables. @@ -82,7 +82,7 @@ class log_Gaussian(prior): return np.exp(np.random.randn(n) * self.sigma + self.mu) -class multivariate_Gaussian: +class MultivariateGaussian: """ Implementation of the multivariate Gaussian probability function, coupled with random variables. @@ -133,20 +133,10 @@ class multivariate_Gaussian: def gamma_from_EV(E, V): - """ - Creates an instance of a gamma prior by specifying the Expected value(s) - and Variance(s) of the distribution. - - :param E: expected value - :param V: variance - - """ warnings.warn("use Gamma.from_EV to create Gamma Prior", FutureWarning) - a = np.square(E) / V - b = E / V - return gamma(a, b) + return Gamma.from_EV(E, V) -class gamma(prior): +class Gamma(Prior): """ Implementation of the Gamma probability function, coupled with random variables. @@ -184,8 +174,20 @@ class gamma(prior): def rvs(self, n): return np.random.gamma(scale=1. / self.b, shape=self.a, size=n) + @staticmethod + def from_EV(E, V): + """ + Creates an instance of a Gamma Prior by specifying the Expected value(s) + and Variance(s) of the distribution. + + :param E: expected value + :param V: variance + """ + a = np.square(E) / V + b = E / V + return Gamma(a, b) -class inverse_gamma(prior): +class inverse_gamma(Prior): """ Implementation of the inverse-Gamma probability function, coupled with random variables. diff --git a/GPy/inference/SCG.py b/GPy/inference/SCG.py index 4318c197..1cec0baa 100644 --- a/GPy/inference/SCG.py +++ b/GPy/inference/SCG.py @@ -156,8 +156,8 @@ def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=500, display=True, xto # beta = 1. # TODO: betareset!! nsuccess = 0 elif success: - gamma = np.dot(gradold - gradnew, gradnew) / (mu) - d = gamma * d - gradnew + Gamma = np.dot(gradold - gradnew, gradnew) / (mu) + d = Gamma * d - gradnew else: # If we get here, then we haven't terminated in the given number of # iterations. diff --git a/GPy/likelihoods/EP.py b/GPy/likelihoods/EP.py index 9f410bc1..caa10987 100644 --- a/GPy/likelihoods/EP.py +++ b/GPy/likelihoods/EP.py @@ -165,7 +165,7 @@ class EP(likelihood): """ Posterior approximation: q(f|y) = N(f| mu, Sigma) Sigma = Diag + P*R.T*R*P.T + K - mu = w + P*gamma + mu = w + P*Gamma """ mu = np.zeros(self.N) LLT = Kmm.copy() @@ -255,10 +255,10 @@ class EP(likelihood): """ Posterior approximation: q(f|y) = N(f| mu, Sigma) Sigma = Diag + P*R.T*R*P.T + K - mu = w + P*gamma + mu = w + P*Gamma """ self.w = np.zeros(self.N) - self.gamma = np.zeros(M) + self.Gamma = np.zeros(M) mu = np.zeros(self.N) P = P0.copy() R = R0.copy() @@ -311,10 +311,10 @@ class EP(likelihood): RTR = np.dot(R.T,np.dot(np.eye(M) - Delta_tau/(1.+Delta_tau*Sigma_diag[i]) * np.dot(Rp_i,Rp_i.T),R)) R = jitchol(RTR).T self.w[i] += (Delta_v - Delta_tau*self.w[i])*dii/dtd1 - self.gamma += (Delta_v - Delta_tau*mu[i])*np.dot(RTR,P[i,:].T) + self.Gamma += (Delta_v - Delta_tau*mu[i])*np.dot(RTR,P[i,:].T) RPT = np.dot(R,P.T) Sigma_diag = Diag + np.sum(RPT.T*RPT.T,-1) - mu = self.w + np.dot(P,self.gamma) + mu = self.w + np.dot(P,self.Gamma) self.iterations += 1 #Sigma recomptutation with Cholesky decompositon Iplus_Dprod_i = 1./(1.+ Diag0 * self.tau_tilde) @@ -326,8 +326,8 @@ class EP(likelihood): RPT = np.dot(R,P.T) Sigma_diag = Diag + np.sum(RPT.T*RPT.T,-1) self.w = Diag * self.v_tilde - self.gamma = np.dot(R.T, np.dot(RPT,self.v_tilde)) - mu = self.w + np.dot(P,self.gamma) + self.Gamma = np.dot(R.T, np.dot(RPT,self.v_tilde)) + mu = self.w + np.dot(P,self.Gamma) epsilon_np1 = sum((self.tau_tilde-self.np1[-1])**2)/self.N epsilon_np2 = sum((self.v_tilde-self.np2[-1])**2)/self.N self.np1.append(self.tau_tilde.copy()) diff --git a/GPy/models/FITC.py b/GPy/models/FITC.py index da2c4d84..e8078780 100644 --- a/GPy/models/FITC.py +++ b/GPy/models/FITC.py @@ -197,8 +197,8 @@ class FITC(sparse_GP): self.RPT = np.dot(self.R,self.P.T) self.Sigma = np.diag(self.Diag) + np.dot(self.RPT.T,self.RPT) self.w = self.Diag * self.likelihood.v_tilde - self.gamma = np.dot(self.R.T, np.dot(self.RPT,self.likelihood.v_tilde)) - self.mu = self.w + np.dot(self.P,self.gamma) + self.Gamma = np.dot(self.R.T, np.dot(self.RPT,self.likelihood.v_tilde)) + self.mu = self.w + np.dot(self.P,self.Gamma) """ Make a prediction for the generalized FITC model diff --git a/GPy/models/generalized_FITC.py b/GPy/models/generalized_FITC.py index 9512dceb..145b4b11 100644 --- a/GPy/models/generalized_FITC.py +++ b/GPy/models/generalized_FITC.py @@ -99,8 +99,8 @@ class generalized_FITC(sparse_GP): self.RPT = np.dot(self.R,self.P.T) self.Sigma = np.diag(self.Diag) + np.dot(self.RPT.T,self.RPT) self.w = self.Diag * self.likelihood.v_tilde - self.gamma = np.dot(self.R.T, np.dot(self.RPT,self.likelihood.v_tilde)) - self.mu = self.w + np.dot(self.P,self.gamma) + self.Gamma = np.dot(self.R.T, np.dot(self.RPT,self.likelihood.v_tilde)) + self.mu = self.w + np.dot(self.P,self.Gamma) # Remove extra term from dL_dpsi1 self.dL_dpsi1 -= mdot(self.Lmi.T,Lmipsi1*self.likelihood.precision.flatten().reshape(1,self.N)) diff --git a/GPy/testing/prior_tests.py b/GPy/testing/prior_tests.py index b3e3a440..d3269560 100644 --- a/GPy/testing/prior_tests.py +++ b/GPy/testing/prior_tests.py @@ -15,12 +15,12 @@ class PriorTests(unittest.TestCase): X, y = X[:, None], y[:, None] m = GPy.models.GP_regression(X, y) m.ensure_default_constraints() - lognormal = GPy.priors.log_Gaussian(1, 2) + lognormal = GPy.priors.LogGaussian(1, 2) m.set_prior('rbf', lognormal) m.randomize() self.assertTrue(m.checkgrad()) - def test_gamma(self): + def test_Gamma(self): xmin, xmax = 1, 2.5*np.pi b, C, SNR = 1, 0, 0.1 X = np.linspace(xmin, xmax, 500) @@ -29,8 +29,8 @@ class PriorTests(unittest.TestCase): X, y = X[:, None], y[:, None] m = GPy.models.GP_regression(X, y) m.ensure_default_constraints() - gamma = GPy.priors.gamma(1, 1) - m.set_prior('rbf', gamma) + Gamma = GPy.priors.Gamma(1, 1) + m.set_prior('rbf', Gamma) m.randomize() self.assertTrue(m.checkgrad())