mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-05-10 12:32:40 +02:00
added domains to priors
This commit is contained in:
Max Zwiessele
parent
dbc4bc3f3c
commit
e02804a671
2 changed files with 83 additions and 80 deletions
|
|
@ -6,19 +6,20 @@ import numpy as np
|
||||||
import pylab as pb
|
import pylab as pb
|
||||||
from scipy.special import gammaln, digamma
|
from scipy.special import gammaln, digamma
|
||||||
from ..util.linalg import pdinv
|
from ..util.linalg import pdinv
|
||||||
from GPy.core.domains import UNDEFINED
|
from GPy.core.domains import REAL, POSITIVE
|
||||||
|
import warnings
|
||||||
|
|
||||||
class prior:
|
class prior:
|
||||||
domain = UNDEFINED
|
domain = None
|
||||||
def pdf(self,x):
|
def pdf(self, x):
|
||||||
return np.exp(self.lnpdf(x))
|
return np.exp(self.lnpdf(x))
|
||||||
|
|
||||||
def plot(self):
|
def plot(self):
|
||||||
rvs = self.rvs(1000)
|
rvs = self.rvs(1000)
|
||||||
pb.hist(rvs,100,normed=True)
|
pb.hist(rvs, 100, normed=True)
|
||||||
xmin,xmax = pb.xlim()
|
xmin, xmax = pb.xlim()
|
||||||
xx = np.linspace(xmin,xmax,1000)
|
xx = np.linspace(xmin, xmax, 1000)
|
||||||
pb.plot(xx,self.pdf(xx),'r',linewidth=2)
|
pb.plot(xx, self.pdf(xx), 'r', linewidth=2)
|
||||||
|
|
||||||
|
|
||||||
class Gaussian(prior):
|
class Gaussian(prior):
|
||||||
|
|
@ -31,24 +32,24 @@ class Gaussian(prior):
|
||||||
.. Note:: Bishop 2006 notation is used throughout the code
|
.. Note:: Bishop 2006 notation is used throughout the code
|
||||||
|
|
||||||
"""
|
"""
|
||||||
|
domain = REAL
|
||||||
def __init__(self,mu,sigma):
|
def __init__(self, mu, sigma):
|
||||||
self.mu = float(mu)
|
self.mu = float(mu)
|
||||||
self.sigma = float(sigma)
|
self.sigma = float(sigma)
|
||||||
self.sigma2 = np.square(self.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):
|
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):
|
def lnpdf(self, x):
|
||||||
return self.constant - 0.5*np.square(x-self.mu)/self.sigma2
|
return self.constant - 0.5 * np.square(x - self.mu) / self.sigma2
|
||||||
|
|
||||||
def lnpdf_grad(self,x):
|
def lnpdf_grad(self, x):
|
||||||
return -(x-self.mu)/self.sigma2
|
return -(x - self.mu) / self.sigma2
|
||||||
|
|
||||||
def rvs(self,n):
|
def rvs(self, n):
|
||||||
return np.random.randn(n)*self.sigma + self.mu
|
return np.random.randn(n) * self.sigma + self.mu
|
||||||
|
|
||||||
|
|
||||||
class log_Gaussian(prior):
|
class log_Gaussian(prior):
|
||||||
|
|
@ -61,24 +62,24 @@ class log_Gaussian(prior):
|
||||||
.. Note:: Bishop 2006 notation is used throughout the code
|
.. Note:: Bishop 2006 notation is used throughout the code
|
||||||
|
|
||||||
"""
|
"""
|
||||||
|
domain = POSITIVE
|
||||||
def __init__(self,mu,sigma):
|
def __init__(self, mu, sigma):
|
||||||
self.mu = float(mu)
|
self.mu = float(mu)
|
||||||
self.sigma = float(sigma)
|
self.sigma = float(sigma)
|
||||||
self.sigma2 = np.square(self.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):
|
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):
|
def lnpdf(self, x):
|
||||||
return self.constant - 0.5*np.square(np.log(x)-self.mu)/self.sigma2 -np.log(x)
|
return self.constant - 0.5 * np.square(np.log(x) - self.mu) / self.sigma2 - np.log(x)
|
||||||
|
|
||||||
def lnpdf_grad(self,x):
|
def lnpdf_grad(self, x):
|
||||||
return -((np.log(x)-self.mu)/self.sigma2+1.)/x
|
return -((np.log(x) - self.mu) / self.sigma2 + 1.) / x
|
||||||
|
|
||||||
def rvs(self,n):
|
def rvs(self, n):
|
||||||
return np.exp(np.random.randn(n)*self.sigma + self.mu)
|
return np.exp(np.random.randn(n) * self.sigma + self.mu)
|
||||||
|
|
||||||
|
|
||||||
class multivariate_Gaussian:
|
class multivariate_Gaussian:
|
||||||
|
|
@ -91,47 +92,47 @@ class multivariate_Gaussian:
|
||||||
.. Note:: Bishop 2006 notation is used throughout the code
|
.. Note:: Bishop 2006 notation is used throughout the code
|
||||||
|
|
||||||
"""
|
"""
|
||||||
|
domain = REAL
|
||||||
def __init__(self,mu,var):
|
def __init__(self, mu, var):
|
||||||
self.mu = np.array(mu).flatten()
|
self.mu = np.array(mu).flatten()
|
||||||
self.var = np.array(var)
|
self.var = np.array(var)
|
||||||
assert len(self.var.shape)==2
|
assert len(self.var.shape) == 2
|
||||||
assert self.var.shape[0]==self.var.shape[1]
|
assert self.var.shape[0] == self.var.shape[1]
|
||||||
assert self.var.shape[0]==self.mu.size
|
assert self.var.shape[0] == self.mu.size
|
||||||
self.D = self.mu.size
|
self.D = self.mu.size
|
||||||
self.inv, self.hld = pdinv(self.var)
|
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):
|
def summary(self):
|
||||||
raise NotImplementedError
|
raise NotImplementedError
|
||||||
|
|
||||||
def pdf(self,x):
|
def pdf(self, x):
|
||||||
return np.exp(self.lnpdf(x))
|
return np.exp(self.lnpdf(x))
|
||||||
|
|
||||||
def lnpdf(self,x):
|
def lnpdf(self, x):
|
||||||
d = x-self.mu
|
d = x - self.mu
|
||||||
return self.constant - 0.5*np.sum(d*np.dot(d,self.inv),1)
|
return self.constant - 0.5 * np.sum(d * np.dot(d, self.inv), 1)
|
||||||
|
|
||||||
def lnpdf_grad(self,x):
|
def lnpdf_grad(self, x):
|
||||||
d = x-self.mu
|
d = x - self.mu
|
||||||
return -np.dot(self.inv,d)
|
return -np.dot(self.inv, d)
|
||||||
|
|
||||||
def rvs(self,n):
|
def rvs(self, n):
|
||||||
return np.random.multivariate_normal(self.mu, self.var,n)
|
return np.random.multivariate_normal(self.mu, self.var, n)
|
||||||
|
|
||||||
def plot(self):
|
def plot(self):
|
||||||
if self.D==2:
|
if self.D == 2:
|
||||||
rvs = self.rvs(200)
|
rvs = self.rvs(200)
|
||||||
pb.plot(rvs[:,0],rvs[:,1], 'kx', mew=1.5)
|
pb.plot(rvs[:, 0], rvs[:, 1], 'kx', mew=1.5)
|
||||||
xmin,xmax = pb.xlim()
|
xmin, xmax = pb.xlim()
|
||||||
ymin,ymax = pb.ylim()
|
ymin, ymax = pb.ylim()
|
||||||
xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
|
xx, yy = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
|
||||||
xflat = np.vstack((xx.flatten(),yy.flatten())).T
|
xflat = np.vstack((xx.flatten(), yy.flatten())).T
|
||||||
zz = self.pdf(xflat).reshape(100,100)
|
zz = self.pdf(xflat).reshape(100, 100)
|
||||||
pb.contour(xx,yy,zz,linewidths=2)
|
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)
|
Creates an instance of a gamma prior by specifying the Expected value(s)
|
||||||
and Variance(s) of the distribution.
|
and Variance(s) of the distribution.
|
||||||
|
|
@ -140,10 +141,10 @@ def gamma_from_EV(E,V):
|
||||||
:param V: variance
|
:param V: variance
|
||||||
|
|
||||||
"""
|
"""
|
||||||
|
warnings.warn("use Gamma.from_EV to create Gamma Prior", FutureWarning)
|
||||||
a = np.square(E)/V
|
a = np.square(E) / V
|
||||||
b = E/V
|
b = E / V
|
||||||
return gamma(a,b)
|
return gamma(a, b)
|
||||||
|
|
||||||
class gamma(prior):
|
class gamma(prior):
|
||||||
"""
|
"""
|
||||||
|
|
@ -155,33 +156,34 @@ class gamma(prior):
|
||||||
.. Note:: Bishop 2006 notation is used throughout the code
|
.. 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.a = float(a)
|
||||||
self.b = float(b)
|
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):
|
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):
|
def summary(self):
|
||||||
ret = {"E[x]": self.a/self.b,\
|
ret = {"E[x]": self.a / self.b, \
|
||||||
"E[ln x]": digamma(self.a) - np.log(self.b),\
|
"E[ln x]": digamma(self.a) - np.log(self.b), \
|
||||||
"var[x]": self.a/self.b/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}
|
"Entropy": gammaln(self.a) - (self.a - 1.) * digamma(self.a) - np.log(self.b) + self.a}
|
||||||
if self.a >1:
|
if self.a > 1:
|
||||||
ret['Mode'] = (self.a-1.)/self.b
|
ret['Mode'] = (self.a - 1.) / self.b
|
||||||
else:
|
else:
|
||||||
ret['mode'] = np.nan
|
ret['mode'] = np.nan
|
||||||
return ret
|
return ret
|
||||||
|
|
||||||
def lnpdf(self,x):
|
def lnpdf(self, x):
|
||||||
return self.constant + (self.a-1)*np.log(x) - self.b*x
|
return self.constant + (self.a - 1) * np.log(x) - self.b * x
|
||||||
|
|
||||||
def lnpdf_grad(self,x):
|
def lnpdf_grad(self, x):
|
||||||
return (self.a-1.)/x - self.b
|
return (self.a - 1.) / x - self.b
|
||||||
|
|
||||||
def rvs(self,n):
|
def rvs(self, n):
|
||||||
return np.random.gamma(scale=1./self.b,shape=self.a,size=n)
|
return np.random.gamma(scale=1. / self.b, shape=self.a, size=n)
|
||||||
|
|
||||||
class inverse_gamma(prior):
|
class inverse_gamma(prior):
|
||||||
"""
|
"""
|
||||||
|
|
@ -193,19 +195,20 @@ class inverse_gamma(prior):
|
||||||
.. Note:: Bishop 2006 notation is used throughout the code
|
.. 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.a = float(a)
|
||||||
self.b = float(b)
|
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):
|
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):
|
def lnpdf(self, x):
|
||||||
return self.constant - (self.a+1)*np.log(x) - self.b/x
|
return self.constant - (self.a + 1) * np.log(x) - self.b / x
|
||||||
|
|
||||||
def lnpdf_grad(self,x):
|
def lnpdf_grad(self, x):
|
||||||
return -(self.a+1.)/x + self.b/x**2
|
return -(self.a + 1.) / x + self.b / x ** 2
|
||||||
|
|
||||||
def rvs(self,n):
|
def rvs(self, n):
|
||||||
return 1./np.random.gamma(scale=1./self.b,shape=self.a,size=n)
|
return 1. / np.random.gamma(scale=1. / self.b, shape=self.a, size=n)
|
||||||
|
|
|
||||||
|
|
@ -3,10 +3,10 @@
|
||||||
|
|
||||||
|
|
||||||
import numpy as np
|
import numpy as np
|
||||||
from GPy.core.domains import UNDEFINED, POSITIVE, NEGATIVE, BOUNDED
|
from GPy.core.domains import POSITIVE, NEGATIVE, BOUNDED
|
||||||
|
|
||||||
class transformation(object):
|
class transformation(object):
|
||||||
domain = UNDEFINED
|
domain = None
|
||||||
def f(self, x):
|
def f(self, x):
|
||||||
raise NotImplementedError
|
raise NotImplementedError
|
||||||
|
|
||||||
|
|
|
||||||
Loading…
Add table
Add a link
Reference in a new issue