From e02804a6712b891d3de67af2ae4f206ddf13b729 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Tue, 4 Jun 2013 17:03:29 +0100 Subject: [PATCH] added domains to priors --- GPy/core/priors.py | 159 ++++++++++++++++++------------------ GPy/core/transformations.py | 4 +- 2 files changed, 83 insertions(+), 80 deletions(-) diff --git a/GPy/core/priors.py b/GPy/core/priors.py index 33bcdc69..74ca63bf 100644 --- a/GPy/core/priors.py +++ b/GPy/core/priors.py @@ -6,19 +6,20 @@ import numpy as np import pylab as pb from scipy.special import gammaln, digamma from ..util.linalg import pdinv -from GPy.core.domains import UNDEFINED +from GPy.core.domains import REAL, POSITIVE +import warnings class prior: - domain = UNDEFINED - def pdf(self,x): + domain = None + def pdf(self, x): return np.exp(self.lnpdf(x)) def plot(self): rvs = self.rvs(1000) - pb.hist(rvs,100,normed=True) - xmin,xmax = pb.xlim() - xx = np.linspace(xmin,xmax,1000) - pb.plot(xx,self.pdf(xx),'r',linewidth=2) + pb.hist(rvs, 100, normed=True) + xmin, xmax = pb.xlim() + xx = np.linspace(xmin, xmax, 1000) + pb.plot(xx, self.pdf(xx), 'r', linewidth=2) class Gaussian(prior): @@ -31,24 +32,24 @@ class Gaussian(prior): .. Note:: Bishop 2006 notation is used throughout the code """ - - def __init__(self,mu,sigma): + domain = REAL + def __init__(self, mu, sigma): self.mu = float(mu) self.sigma = float(sigma) self.sigma2 = np.square(self.sigma) - self.constant = -0.5*np.log(2*np.pi*self.sigma2) + self.constant = -0.5 * np.log(2 * np.pi * self.sigma2) def __str__(self): - return "N("+str(np.round(self.mu))+', '+str(np.round(self.sigma2))+')' + return "N(" + str(np.round(self.mu)) + ', ' + str(np.round(self.sigma2)) + ')' - def lnpdf(self,x): - return self.constant - 0.5*np.square(x-self.mu)/self.sigma2 + def lnpdf(self, x): + return self.constant - 0.5 * np.square(x - self.mu) / self.sigma2 - def lnpdf_grad(self,x): - return -(x-self.mu)/self.sigma2 + def lnpdf_grad(self, x): + return -(x - self.mu) / self.sigma2 - def rvs(self,n): - return np.random.randn(n)*self.sigma + self.mu + def rvs(self, n): + return np.random.randn(n) * self.sigma + self.mu class log_Gaussian(prior): @@ -61,24 +62,24 @@ class log_Gaussian(prior): .. Note:: Bishop 2006 notation is used throughout the code """ - - def __init__(self,mu,sigma): + domain = POSITIVE + def __init__(self, mu, sigma): self.mu = float(mu) self.sigma = float(sigma) self.sigma2 = np.square(self.sigma) - self.constant = -0.5*np.log(2*np.pi*self.sigma2) + self.constant = -0.5 * np.log(2 * np.pi * self.sigma2) def __str__(self): - return "lnN("+str(np.round(self.mu))+', '+str(np.round(self.sigma2))+')' + return "lnN(" + str(np.round(self.mu)) + ', ' + str(np.round(self.sigma2)) + ')' - def lnpdf(self,x): - return self.constant - 0.5*np.square(np.log(x)-self.mu)/self.sigma2 -np.log(x) + def lnpdf(self, x): + return self.constant - 0.5 * np.square(np.log(x) - self.mu) / self.sigma2 - np.log(x) - def lnpdf_grad(self,x): - return -((np.log(x)-self.mu)/self.sigma2+1.)/x + def lnpdf_grad(self, x): + return -((np.log(x) - self.mu) / self.sigma2 + 1.) / x - def rvs(self,n): - return np.exp(np.random.randn(n)*self.sigma + self.mu) + def rvs(self, n): + return np.exp(np.random.randn(n) * self.sigma + self.mu) class multivariate_Gaussian: @@ -91,47 +92,47 @@ class multivariate_Gaussian: .. Note:: Bishop 2006 notation is used throughout the code """ - - def __init__(self,mu,var): + domain = REAL + def __init__(self, mu, var): self.mu = np.array(mu).flatten() self.var = np.array(var) - assert len(self.var.shape)==2 - assert self.var.shape[0]==self.var.shape[1] - assert self.var.shape[0]==self.mu.size + assert len(self.var.shape) == 2 + assert self.var.shape[0] == self.var.shape[1] + assert self.var.shape[0] == self.mu.size self.D = self.mu.size self.inv, self.hld = pdinv(self.var) - self.constant = -0.5*self.D*np.log(2*np.pi) - self.hld + self.constant = -0.5 * self.D * np.log(2 * np.pi) - self.hld def summary(self): raise NotImplementedError - def pdf(self,x): + def pdf(self, x): return np.exp(self.lnpdf(x)) - def lnpdf(self,x): - d = x-self.mu - return self.constant - 0.5*np.sum(d*np.dot(d,self.inv),1) + def lnpdf(self, x): + d = x - self.mu + return self.constant - 0.5 * np.sum(d * np.dot(d, self.inv), 1) - def lnpdf_grad(self,x): - d = x-self.mu - return -np.dot(self.inv,d) + def lnpdf_grad(self, x): + d = x - self.mu + return -np.dot(self.inv, d) - def rvs(self,n): - return np.random.multivariate_normal(self.mu, self.var,n) + def rvs(self, n): + return np.random.multivariate_normal(self.mu, self.var, n) def plot(self): - if self.D==2: + if self.D == 2: rvs = self.rvs(200) - pb.plot(rvs[:,0],rvs[:,1], 'kx', mew=1.5) - xmin,xmax = pb.xlim() - ymin,ymax = pb.ylim() + pb.plot(rvs[:, 0], rvs[:, 1], 'kx', mew=1.5) + xmin, xmax = pb.xlim() + ymin, ymax = pb.ylim() xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j] - xflat = np.vstack((xx.flatten(),yy.flatten())).T - zz = self.pdf(xflat).reshape(100,100) - pb.contour(xx,yy,zz,linewidths=2) + xflat = np.vstack((xx.flatten(), yy.flatten())).T + zz = self.pdf(xflat).reshape(100, 100) + pb.contour(xx, yy, zz, linewidths=2) -def gamma_from_EV(E,V): +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. @@ -140,10 +141,10 @@ def gamma_from_EV(E,V): :param V: variance """ - - a = np.square(E)/V - b = E/V - return gamma(a,b) + warnings.warn("use Gamma.from_EV to create Gamma Prior", FutureWarning) + a = np.square(E) / V + b = E / V + return gamma(a, b) class gamma(prior): """ @@ -155,33 +156,34 @@ class gamma(prior): .. Note:: Bishop 2006 notation is used throughout the code """ - def __init__(self,a,b): + domain = POSITIVE + def __init__(self, a, b): self.a = float(a) self.b = float(b) - self.constant = -gammaln(self.a) + a*np.log(b) + self.constant = -gammaln(self.a) + a * np.log(b) def __str__(self): - return "Ga("+str(np.round(self.a))+', '+str(np.round(self.b))+')' + return "Ga(" + str(np.round(self.a)) + ', ' + str(np.round(self.b)) + ')' def summary(self): - ret = {"E[x]": self.a/self.b,\ - "E[ln x]": digamma(self.a) - np.log(self.b),\ - "var[x]": self.a/self.b/self.b,\ - "Entropy": gammaln(self.a) - (self.a-1.)*digamma(self.a) - np.log(self.b) + self.a} - if self.a >1: - ret['Mode'] = (self.a-1.)/self.b + ret = {"E[x]": self.a / self.b, \ + "E[ln x]": digamma(self.a) - np.log(self.b), \ + "var[x]": self.a / self.b / self.b, \ + "Entropy": gammaln(self.a) - (self.a - 1.) * digamma(self.a) - np.log(self.b) + self.a} + if self.a > 1: + ret['Mode'] = (self.a - 1.) / self.b else: ret['mode'] = np.nan return ret - def lnpdf(self,x): - return self.constant + (self.a-1)*np.log(x) - self.b*x + def lnpdf(self, x): + return self.constant + (self.a - 1) * np.log(x) - self.b * x - def lnpdf_grad(self,x): - return (self.a-1.)/x - self.b + def lnpdf_grad(self, x): + return (self.a - 1.) / x - self.b - def rvs(self,n): - return np.random.gamma(scale=1./self.b,shape=self.a,size=n) + def rvs(self, n): + return np.random.gamma(scale=1. / self.b, shape=self.a, size=n) class inverse_gamma(prior): """ @@ -193,19 +195,20 @@ class inverse_gamma(prior): .. Note:: Bishop 2006 notation is used throughout the code """ - def __init__(self,a,b): + domain = POSITIVE + def __init__(self, a, b): self.a = float(a) self.b = float(b) - self.constant = -gammaln(self.a) + a*np.log(b) + self.constant = -gammaln(self.a) + a * np.log(b) def __str__(self): - return "iGa("+str(np.round(self.a))+', '+str(np.round(self.b))+')' + return "iGa(" + str(np.round(self.a)) + ', ' + str(np.round(self.b)) + ')' - def lnpdf(self,x): - return self.constant - (self.a+1)*np.log(x) - self.b/x + def lnpdf(self, x): + return self.constant - (self.a + 1) * np.log(x) - self.b / x - def lnpdf_grad(self,x): - return -(self.a+1.)/x + self.b/x**2 + def lnpdf_grad(self, x): + return -(self.a + 1.) / x + self.b / x ** 2 - def rvs(self,n): - return 1./np.random.gamma(scale=1./self.b,shape=self.a,size=n) + def rvs(self, n): + return 1. / np.random.gamma(scale=1. / self.b, shape=self.a, size=n) diff --git a/GPy/core/transformations.py b/GPy/core/transformations.py index b9748984..2520a33b 100644 --- a/GPy/core/transformations.py +++ b/GPy/core/transformations.py @@ -3,10 +3,10 @@ import numpy as np -from GPy.core.domains import UNDEFINED, POSITIVE, NEGATIVE, BOUNDED +from GPy.core.domains import POSITIVE, NEGATIVE, BOUNDED class transformation(object): - domain = UNDEFINED + domain = None def f(self, x): raise NotImplementedError