some love for the priors class

This commit is contained in:
Nicolò Fusi 2013-01-27 18:23:19 +00:00
parent c43904c3bf
commit b73cb33a94

View file

@ -10,6 +10,7 @@ from ..util.linalg import pdinv
class prior:
def pdf(self,x):
return np.exp(self.lnpdf(x))
def plot(self):
rvs = self.rvs(1000)
pb.hist(rvs,100,normed=True)
@ -17,45 +18,79 @@ class prior:
xx = np.linspace(xmin,xmax,1000)
pb.plot(xx,self.pdf(xx),'r',linewidth=2)
class Gaussian(prior):
"""
Implementation of the univariate Gaussian probability function, coupled with random variables, since scipy.stats sucks.
Using Bishop 2006 notation"""
Implementation of the univariate Gaussian probability function, coupled with random variables.
:param mu: mean
:param sigma: standard deviation
.. Note:: Bishop 2006 notation is used throughout the code
"""
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)
def __str__(self):
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_grad(self,x):
return -(x-self.mu)/self.sigma2
def rvs(self,n):
return np.random.randn(n)*self.sigma + self.mu
class log_Gaussian(prior):
"""
Implementation of the univariate *log*-Gaussian probability function, coupled with random variables.
:param mu: mean
:param sigma: standard deviation
.. Note:: Bishop 2006 notation is used throughout the code
"""
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)
def __str__(self):
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_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)
class multivariate_Gaussian:
"""
Implementation of the multivariate Gaussian probability function, coupled with random variables, since scipy.stats sucks.
Using Bishop 2006 notation"""
Implementation of the multivariate Gaussian probability function, coupled with random variables.
:param mu: mean (N-dimensional array)
:param var: covariance matrix (NxN)
.. Note:: Bishop 2006 notation is used throughout the code
"""
def __init__(self,mu,var):
self.mu = np.array(mu).flatten()
self.var = np.array(var)
@ -67,17 +102,22 @@ class multivariate_Gaussian:
self.constant = -0.5*self.D*np.log(2*np.pi) - self.hld
def summary(self):
pass #TODO
raise NotImplementedError
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_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 plot(self):
if self.D==2:
rvs = self.rvs(200)
@ -91,7 +131,15 @@ class multivariate_Gaussian:
def gamma_from_EV(E,V):
"""create an instance of a gamma prior by specifying the Expected value(s) and Variance(s) of the distribution"""
"""
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)
@ -99,15 +147,23 @@ def gamma_from_EV(E,V):
class gamma(prior):
"""
Implementation of the Gamma probability function, coupled with random variables, since scipy.stats sucks.
Using Bishop 2006 notation
Implementation of the Gamma probability function, coupled with random variables.
:param a: shape parameter
:param b: rate parameter (warning: it's the *inverse* of the scale)
.. Note:: Bishop 2006 notation is used throughout the code
"""
def __init__(self,a,b):
self.a = float(a)
self.b = float(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))+')'
def summary(self):
ret = {"E[x]": self.a/self.b,\
"E[ln x]": digamma(self.a) - np.log(self.b),\
@ -118,10 +174,12 @@ class gamma(prior):
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_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)