added calc_utils

Author: jsalvatier

PyMC.tmbundle/Snippets/Impute.tmSnippet

@@ -0,0 +1,31 @@
+'''
+Created on Jan 20, 2011
+
+@author: jsalvatier
+'''
+import numpy as np
+from collections import defaultdict
+
+_sts_memory = defaultdict(dict)
+def sum_to_shape(key1,key2, value, sum_shape):
+    
+    try :
+        axes, lx = _sts_memory[key1][key2]
+    
+    except KeyError: 
+        
+        value_shape = np.array(np.shape(value))
+        
+        sum_shape_expanded = np.zeros(value_shape.size)
+        sum_shape_expanded[0:len(sum_shape)] += np.array(sum_shape)
+        
+        axes = np.where(sum_shape_expanded != value_shape)[0]
+        lx = np.size(axes)
+        
+        _sts_memory[key1][key2] = (axes, lx )
+         
+    if lx > 0:
+        return np.apply_over_axes(np.sum, value, axes)
+
+    else:
+        return value

pymc/calc_utils.py

@@ -59,8 +59,8 @@ class ArgumentError(AttributeError):
                                'degenerate', 'exponential', 'exponweib',
                                'gamma', 'half_normal', 'hypergeometric',
                                'inverse_gamma', 'laplace', 'logistic',
-                               'lognormal', 'normal', 'pareto', 't', 
-                               'truncated_pareto', 'uniform',
+                               'lognormal', 'noncentral_t', 'normal', 
+                               'pareto', 't', 'truncated_pareto', 'uniform',
                                'weibull', 'skew_normal', 'truncated_normal',
                                'von_mises']
 sc_bool_distributions = ['bernoulli']
@@ -2241,9 +2241,9 @@ def truncated_pareto_expval(alpha, m, b):
 
     if alpha <= 1:
         return inf
-    part1 = (m**alpha)/(1 - (m/b)**alpha)
-    part2 = alpha/(alpha-1)
-    part3 = (1/(m**(alpha-1)) - 1/(b**(alpha-1)))
+    part1 = (m**alpha)/(1. - (m/b)**alpha)
+    part2 = 1.*alpha/(alpha-1)
+    part3 = (1./(m**(alpha-1)) - 1./(b**(alpha-1.)))
     return part1*part2*part3
 
 def truncated_pareto_like(x, alpha, m, b):
@@ -2404,72 +2404,23 @@ def truncated_poisson_like(x,mu,k):
 
 # Truncated normal distribution--------------------------
 @randomwrap
-def rtruncated_normal(mu, tau, a=None, b=None, size=None):
+def rtruncated_normal(mu, tau, a=-np.inf, b=np.inf, size=None):
     """rtruncated_normal(mu, tau, a, b, size=1)
 
-    Random truncated normal variates using method from Robert (1995).
+    Random truncated normal variates.
     """
-
-    factor = 10
-    sign = 1.0
-
-    while True:
-
-        if a is None and b is None:
-            raise ValueError, 'No truncation boundary given.'
-
-        elif a is None or b is None:
-            # One-sided truncation
-            
-            if a is None:
-                # See top of p.123 in Robert (1995)
-                a = -b
-                mu = -mu
-                sign = -1.0
-            
-            # Algorithm is in terms of standard normal
-            a = np.sqrt(tau) * (a - mu)
-            
-            # Parameter of exponential proposal
-            beta = (a + np.sqrt(a**2 + 4))/2.0
-            # Sample from exponential
-            z = np.random.exponential(1./beta, size*factor) + a
-
-            if a<beta:
-
-                x = np.exp(-0.5 * (beta - z)**2)
-
-            else:
-                
-                x = np.exp(0.5 * (a - beta)**2) * np.exp(-(beta - z)**2)
-
-        else:
-            # Two-sided truncation
-
-            # Algorithm is in terms of standard normal
-            a = np.sqrt(tau) * (a - mu)
-            b = np.sqrt(tau) * (b - mu)
-
-            # Sample z ~ U(a,b)
-            z = (b - a) * random_number(size*factor) + a
-            
-            if a<=0<=b:
-                x = np.exp(-0.5 * z**2)
-            elif b<0:
-                x = np.exp(0.5 * (b**2 - z**2))
-            else:
-                x = np.exp(0.5 * (a**2 - z**2))
-
-        # Accept-reject
-        u = random_number(size*factor)
-        y = sign * (z[u <= x] / np.sqrt(tau) + mu)
-
-        # Return <size> samples
-        if len(y) >= size:
-            return y[:size]
-        else:
-            # Get a larger sample next time
-            factor *=10
+    
+    sigma = 1./np.sqrt(tau)
+    na = pymc.utils.normcdf((a-mu)/sigma)
+    nb = pymc.utils.normcdf((b-mu)/sigma)
+    
+    # Use the inverse CDF generation method.
+    U = np.random.mtrand.uniform(size=size)	
+    q = U * nb + (1-U)*na
+    R = pymc.utils.invcdf(q)
+    
+    # Unnormalize
+    return R*sigma + mu
 
 rtruncnorm = rtruncated_normal
 
@@ -2496,7 +2447,7 @@ def truncated_normal_expval(mu, tau, a, b):
 
 truncnorm_expval = truncated_normal_expval
 
-def truncated_normal_like(x, mu, tau, a, b):
+def truncated_normal_like(x, mu, tau, a=None, b=None):
     R"""truncnorm_like(x, mu, tau, a, b)
 
     Truncated normal log-likelihood.
@@ -2514,7 +2465,9 @@ def truncated_normal_like(x, mu, tau, a, b):
       - `b` : Right bound of the distribution.
     """
     x = np.atleast_1d(x)
+    if a is None: a = -np.inf
     a = np.atleast_1d(a)
+    if b is None: b = np.inf
     b = np.atleast_1d(b)
     mu = np.atleast_1d(mu)
     sigma = (1./np.atleast_1d(np.sqrt(tau)))
@@ -2617,6 +2570,50 @@ def t_expval(nu):
     Expectation of Student's t random variables.
     """
     return 0
+    
+# Non-central Student's t-----------------------------------
+@randomwrap
+def rnoncentral_t(mu, lam, nu, size=None):
+    """rnoncentral_t(mu, lam, nu, size=1)
+
+    Non-central Student's t random variates.
+    """
+    tau = rgamma(nu/2., nu/(2.*lam), size)
+    return rnormal(mu, tau)
+
+def noncentral_t_like(x, mu, lam, nu):
+    R"""noncentral_t_like(x, mu, lam, nu)
+
+    Non-central Student's T log-likelihood. Describes a normal variable
+    whose precision is gamma distributed.
+
+    .. math::
+        f(x|\mu,\lambda,\nu) = \frac{\Gamma(\frac{\nu +
+        1}{2})}{\Gamma(\frac{\nu}{2})}
+        \left(\frac{\lambda}{\pi\nu}\right)^{\frac{1}{2}}
+        \left[1+\frac{\lambda(x-\mu)^2}{\nu}\right]^{-\frac{\nu+1}{2}}
+
+    :Parameters:
+      - `x` : Input data.
+      - `mu` : Location parameter.
+      - `lambda` : Scale parameter. 
+      - `nu` : Degrees of freedom.
+
+    """
+    mu = np.asarray(mu)
+    lam = np.asarray(lam)
+    nu = np.asarray(nu)
+    return flib.nct(x, mu, lam, nu)
+
+def noncentral_t_expval(mu, lam, nu):
+    """noncentral_t_expval(mu, lam, nu)
+
+    Expectation of non-central Student's t random variables. Only defined
+    for nu>1.
+    """
+    if nu>1:
+        return mu
+    return inf
 
 def t_grad_setup(x, nu, f):
     nu = np.asarray(nu)

pymc/distributions.py

@@ -1443,6 +1443,65 @@ SUBROUTINE chi2_grad_nu(x,nu,n,nnu,gradlikenu)
       return
       END
 
+      SUBROUTINE nct(x,mu,lam,nu,n,nmu,nlam,nnu,like)
+
+c Non-central Student's t log-likelihood function    
+
+cf2py double precision dimension(n),intent(in) :: x
+cf2py double precision dimension(nmu),intent(in) :: mu
+cf2py double precision dimension(nlam),intent(in) :: lam
+cf2py double precision dimension(nnu),intent(in) :: nu
+cf2py double precision intent(out) :: like
+cf2py integer intent(hide),depend(x) :: n=len(x)
+cf2py integer intent(hide),depend(mu) :: nmu=len(mu)
+cf2py integer intent(hide),depend(lam) :: nlam=len(lam)
+cf2py integer intent(hide),depend(nu) :: nnu=len(nu)
+cf2py threadsafe
+
+      IMPLICIT NONE
+      INTEGER n, i, nnu, nmu, nlam
+      DOUBLE PRECISION x(n)
+      DOUBLE PRECISION nu(nnu), mu(nmu), lam(nlam), like, infinity
+      DOUBLE PRECISION mut, lamt, nut
+      PARAMETER (infinity = 1.7976931348623157d308)
+      DOUBLE PRECISION gammln
+      DOUBLE PRECISION PI
+      PARAMETER (PI=3.141592653589793238462643d0)
+      
+      nut = nu(1)
+      mut = mu(1)
+      lamt = lam(1)
+
+      like = 0.0
+      do i=1,n
+        if (nmu .GT. 1) then
+          mut = mu(i)
+        endif
+        if (nlam .GT. 1) then
+          lamt = lam(i)
+        endif
+        if (nnu .GT. 1) then
+          nut = nu(i)
+        endif
+        
+        if (nut .LE. 0.0) then
+          like = -infinity
+          RETURN
+        endif
+        if (lamt .LE. 0.0) then
+          like = -infinity
+          RETURN
+        endif
+
+        like = like + gammln((nut+1.0)/2.0)
+        like = like - gammln(nut/2.0)
+        like = like + 0.5*dlog(lamt) - 0.5*dlog(nut * PI)
+        like = like - (nut+1)/2 * dlog(1 + (lamt*(x(i) - mut)**2)/nut)
+      enddo
+      return
+      END
+
+
       SUBROUTINE multinomial(x,n,p,nx,nn,np,k,like)
 
 c Multinomial log-likelihood function     

pymc/flib.f

@@ -511,6 +511,25 @@ def test_consistency(self):
             compare_hist(figname='Student t', **figdata)
         assert_array_almost_equal(hist, like,1)
 
+class test_noncentral_t(TestCase):
+    """Based on gamma."""
+    def test_consistency(self):
+        parameters={'mu':-10, 'lam':0.2, 'nu':5}
+        hist, like, figdata = consistency(rnoncentral_t, noncentral_t_like,
+            parameters, nrandom=5000)
+        if PLOT:
+            compare_hist(figname='noncentral t', **figdata)
+        assert_array_almost_equal(hist, like,1)
+        
+    def test_vectorization(self):
+        a = flib.nct([3,4,5], mu=3, lam=.1, nu=5)
+        b = flib.nct([3,4,5], mu=[3,3,3], lam=.1, nu=5)
+        c = flib.nct([3,4,5], mu=[3,3,3], lam=[.1,.1,.1], nu=5)
+        d = flib.nct([3,4,5], mu=[3,3,3], lam=[.1,.1,.1], nu=[5,5,5])
+        assert_equal(a,b)
+        assert_equal(b,c)
+        assert_equal(c,d)
+
 class test_exponweib(TestCase):
     def test_consistency(self):
         parameters = {'alpha':2, 'k':2, 'loc':1, 'scale':3}
@@ -823,6 +842,12 @@ def test_consistency(self):
         if PLOT:
             compare_hist(figname='truncated_pareto', **figdata)
         assert_array_almost_equal(hist, like, 1)
+        
+    def test_random(self):
+        r = rtruncated_pareto(alpha=3, m=1, b=6, size=10000)
+        assert_almost_equal(r.mean(), truncated_pareto_expval(3, 1, 6), 1)
+        assert (r > 1).all()
+        assert (r < 6).all()
 
     def test_vectorization(self):
         a = flib.truncated_pareto([3,4,5], alpha=3, m=1, b=6)
@@ -856,7 +881,12 @@ def test_consistency(self):
         if PLOT:
             compare_hist(figname='poisson', **figdata)
         assert_array_almost_equal(hist, like,1)
-
+        
+    def test_random(self):
+        r = rtruncated_poisson(mu=5, k=1, size=10000)
+        assert_almost_equal(r.mean(), truncated_poisson_expval(5, 1), 1)
+        assert (r >= 1).all()
+        
     def test_normalization(self):
         parameters = {'mu':4., 'k':1}
         summation=discrete_normalization(flib.trpoisson,parameters,20)