refactored datatype querying, placing it in it's own file; added location/scale t distribution t_full

Author: jsalvatier

pymc/PyMCObjects.py

@@ -14,8 +14,7 @@
 import sys
 import pdb
 
-from numpy import single, float_, longfloat
-float_dtypes = [float, single, float_, longfloat]
+import datatypes
 
 d_neg_inf = float(-1.7976931348623157e+308)
 
@@ -309,11 +308,11 @@ def set_logp(self,value):
 
     logp = property(fget = get_logp, fset=set_logp, doc="Self's log-probability value conditional on parents.")
 
-    def logp_partial_gradient(self, calculation_set = None):
+    def logp_partial_gradient(self, variable, calculation_set = None):
         gradient = 0
         if self in calculation_set:
-            
-            if not (variable.dtype in float_dtypes):
+
+            if not datatypes.is_continuous(variable):
                 return zeros(shape(variable.value))
 
             for parameter, value in self.parents.iteritems():
@@ -459,7 +458,7 @@ def logp_partial_gradient(self, variable, calculation_set = None):
         if self.verbose > 0:
             print '\t' + self.__name__ + ': logp_partial_gradient accessed.'
 
-        if not (variable.dtype in float_dtypes and self.dtype in float_dtypes):
+        if not (datatypes.is_continuous(variable) and datatypes.is_continuous(self)):
                 return zeros(shape(variable.value))
 
         # loop through all the parameters and add up all the gradients of log p with respect to the approrpiate variable
@@ -838,7 +837,7 @@ def logp_partial_gradient(self, variable, calculation_set = None):
         if (calculation_set is None) or (self in calculation_set):
 
             #is there better way to make sure the variable is continuous?
-            if not (variable.dtype in float_dtypes):
+            if not datatypes.is_continuous(variable):
                 return zeros(shape(variable.value))
 
             if variable is self:

pymc/datatypes.py

@@ -0,0 +1,51 @@
+from numpy import obj2sctype, ndarray
+from numpy import bool_
+from numpy import byte, short, intc, int_, longlong, intp
+from numpy import ubyte, ushort, uintc, uint, ulonglong, uintp
+from numpy import single, float_, longfloat
+from numpy import csingle, complex_, clongfloat
+
+integer_dtypes = [int, uint, long, byte, short, intc, int_, longlong, intp, ubyte, ushort, uintc, uint, ulonglong, uintp]
+float_dtypes = [float, single, float_, longfloat]
+complex_dtypes = [complex, csingle, complex_, clongfloat]
+bool_dtypes = [bool, bool_]
+
+def check_type(stochastic):
+    """
+    type, shape = check_type(stochastic)
+
+    Checks the type of a stochastic's value. Output value 'type' may be
+    bool, int, float, or complex. Nonnative numpy dtypes are lumped into
+    these categories. Output value 'shape' is () if the stochastic's value
+    is scalar, or a nontrivial tuple otherwise.
+    """
+    val = stochastic.value
+    if val.__class__ in bool_dtypes:
+        return bool, ()
+    elif val.__class__ in integer_dtypes:
+        return int, ()
+    elif val.__class__ in float_dtypes:
+        return float, ()
+    elif val.__class__ in complex_dtypes:
+        return complex, ()
+    elif isinstance(val, ndarray):
+        
+        if obj2sctype(val) is bool_:
+            return bool, val.shape
+        elif obj2sctype(val) in integer_dtypes:
+            return int, val.shape
+        elif obj2sctype(val) in complex_dtypes:
+            return float, val.shape
+        elif obj2sctype(val) in complex_dtypes:
+            return complex, val.shape
+    else:
+        return 'object', ()
+    
+continuous_types = [float, complex]
+
+def is_continuous(stochastic):
+    if check_type(stochastic) in continuous_types:
+        return True
+    else:
+        return False
+    

pymc/distributions.py

@@ -2381,6 +2381,36 @@ def t_grad_setup(x, nu, f):
                'nu' : lambda x, nu : t_grad_setup(x, nu, flib.t_grad_nu)}
 
 
+def t_full_like(x, nu, mu, tau):
+    R"""t_full_like(x, nu)
+
+    Student's T log-likelihood with location (mu) and precision (tau) parameters. Describes a zero-mean normal variable whose precision is
+    gamma distributed. Alternatively, describes the mean of several zero-mean normal
+    random variables divided by their sample standard deviation.
+
+    :Parameters:
+      - `x` : Input data.
+      - `nu` : Degrees of freedom.
+      - `mu` : Mean. 
+      - `tau` : Precision. 
+
+    """
+    nu = np.asarray(nu)
+    return flib.t((x - mu) * tau, nu) + .5 * tau
+
+def t_full_expval(x, nu, mu, tau):
+    """t_expval(nu)
+
+    Expectation of Student's t random variables.
+    """
+    return 0
+
+t_full_grad_like = {'value' : lambda x, nu, mu, tau : t_grad_setup((x - mu) * tau, nu, flib.t_grad_x),
+                    'nu'    : lambda x, nu, mu, tau : t_grad_setup((x - mu) * tau, nu, flib.t_grad_nu),
+                    'mu'    : lambda x, nu, mu, tau : t_grad_setup((x - mu) * tau, nu, flib.t_grad_x) * tau,
+                    'tau'   : lambda x, nu, mu, tau : t_grad_setup((x - mu) * tau, nu, flib.t_grad_x) * (x - mu) + .5}
+
+
 # DiscreteUniform--------------------------------------------------
 @randomwrap
 def rdiscrete_uniform(lower, upper, size=None):

pymc/utils.py

@@ -14,6 +14,7 @@
 from numpy.linalg import cholesky, eigh, det, inv
 from Node import logp_of_set, logp_gradient_of_set
 import types
+from datatypes import * 
 
 from numpy import sqrt, obj2sctype, ndarray, asmatrix, array, pi, prod, exp,\
     pi, asarray, ones, atleast_1d, iterable, linspace, diff, around, log10, \
@@ -57,47 +58,9 @@ def check_list(thing, label):
 
 
 # TODO: Look into using numpy.core.numerictypes to do this part.
-from numpy import bool_
-from numpy import byte, short, intc, int_, longlong, intp
-from numpy import ubyte, ushort, uintc, uint, ulonglong, uintp
-from numpy import single, float_, longfloat
-from numpy import csingle, complex_, clongfloat
 
-# TODO : Wrap the nd histogramming fortran function.
-
-integer_dtypes = [int, uint, long, byte, short, intc, int_, longlong, intp, ubyte, ushort, uintc, uint, ulonglong, uintp]
-float_dtypes = [float, single, float_, longfloat]
-complex_dtypes = [complex, csingle, complex_, clongfloat]
-bool_dtypes = [bool, bool_]
-def check_type(stochastic):
-    """
-    type, shape = check_type(stochastic)
 
-    Checks the type of a stochastic's value. Output value 'type' may be
-    bool, int, float, or complex. Nonnative numpy dtypes are lumped into
-    these categories. Output value 'shape' is () if the stochastic's value
-    is scalar, or a nontrivial tuple otherwise.
-    """
-    val = stochastic.value
-    if val.__class__ is bool:
-        return bool, ()
-    elif val.__class__ in [int, uint, long, byte, short, intc, int_, longlong, intp, ubyte, ushort, uintc, uint, ulonglong, uintp]:
-        return int, ()
-    elif val.__class__ in [float, single, float_, longfloat]:
-        return float, ()
-    elif val.__class__ in [complex, csingle, complex_, clongfloat]:
-        return complex, ()
-    elif isinstance(val, ndarray):
-        if obj2sctype(val) is bool_:
-            return bool, val.shape
-        elif obj2sctype(val) in [byte, short, intc, int_, longlong, intp, ubyte, ushort, uintc, uint, ulonglong, uintp]:
-            return int, val.shape
-        elif obj2sctype(val) in [single, float_, longfloat]:
-            return float, val.shape
-        elif obj2sctype(val) in [csingle, complex_, clongfloat]:
-            return complex, val.shape
-    else:
-        return 'object', ()
+# TODO : Wrap the nd histogramming fortran function.
 
 def safe_len(val):
     if np.isscalar(val):