make it a package

Author: hildensia

.idea/inspectionProfiles/Project_Default.xml

@@ -0,0 +1,3 @@
+#!/usr/bin/env python2
+
+__version__ = '0.2.dev1'

.idea/inspectionProfiles/profiles_settings.xml

@@ -1,11 +1,18 @@
 from __future__ import division
 import numpy as np
 from scipy.special import gammaln
-from scipy.misc import comb
+from scipy.misc import comb, logsumexp
 from decorator import decorator
 
 
 def _dynamicProgramming(f, *args, **kwargs):
+    if f.data is None:
+        f.data = args[0]
+
+    if not np.array_equal(f.data, args[0]):
+        f.cache = {}
+        f.data = args[0]
+
     try:
         f.cache[args[1:3]]
     except KeyError:
@@ -15,10 +22,13 @@ def _dynamicProgramming(f, *args, **kwargs):
 
 def dynamicProgramming(f):
     f.cache = {}
+    f.data = None
     return decorator(_dynamicProgramming, f)
 
 
-def offline_changepoint_detection(data, prior_func, observation_log_likelihood_function, P=None, truncate=-np.inf):
+def offline_changepoint_detection(data, prior_func,
+                                  observation_log_likelihood_function,
+                                  truncate=-np.inf):
     """Compute the likelihood of changepoints on data.
 
     Keyword arguments:
@@ -34,10 +44,8 @@ def offline_changepoint_detection(data, prior_func, observation_log_likelihood_f
     Q = np.zeros((n,))
     g = np.zeros((n,))
     G = np.zeros((n,))
-    
-    if P is None:
-        P = np.ones((n,n))  # a log_likelihood won't become one
- 
+    P = np.ones((n, n)) * -np.inf
+
     # save everything in log representation
     for t in range(n):
         g[t] = np.log(prior_func(t))
@@ -46,28 +54,24 @@ def offline_changepoint_detection(data, prior_func, observation_log_likelihood_f
         else:
             G[t] = np.logaddexp(G[t-1], g[t])
 
-    if P[n-1, n-1] == 1:
-        P[n-1, n-1] = observation_log_likelihood_function(data, n-1, n)
-        
+    P[n-1, n-1] = observation_log_likelihood_function(data, n-1, n)
     Q[n-1] = P[n-1, n-1]
 
     for t in reversed(range(n-1)):
         P_next_cp = -np.inf  # == -log(0)
         for s in range(t, n-1):
-            # don't recompute log likelihoods already saved
-            if P[t, s] == 1:
-                P[t, s] = observation_log_likelihood_function(data, t, s + 1)
-            
+            P[t, s] = observation_log_likelihood_function(data, t, s + 1)
+
             # compute recursion
             summand = P[t, s] + Q[s + 1] + g[s + 1 - t]
             P_next_cp = np.logaddexp(P_next_cp, summand)
-            
-            # truncate sum to become approx. linear in time (see Fearnhead, 2006, eq. (3))
+
+            # truncate sum to become approx. linear in time (see
+            # Fearnhead, 2006, eq. (3))
             if summand - P_next_cp < truncate:
                 break
-                
-        if P[t, n-1] == 1:
-            P[t, n-1] = observation_log_likelihood_function(data, t, n)
+
+        P[t, n-1] = observation_log_likelihood_function(data, t, n)
 
         # (1 - G) is numerical stable until G becomes numerically 1
         if G[n-1-t] < -1e-15:  # exp(-1e-15) = .99999...
@@ -80,21 +84,17 @@ def offline_changepoint_detection(data, prior_func, observation_log_likelihood_f
 
     Pcp = np.ones((n-1, n-1)) * -np.inf
     for t in range(n-1):
-        if P[0, t] == 1:
-            P[0, t] = observation_log_likelihood_function(data, 0, t+1)
-        Pcp[0, t] = P[0, t] + Q[t + 1] + g[t] - Q[0] 
+        Pcp[0, t] = P[0, t] + Q[t + 1] + g[t] - Q[0]
     for j in range(1, n-1):
         for t in range(j, n-1):
-            for tau in range(j-1, t):
-                if P[tau+1, t] == 1:
-                    P[tau+1, t] = observation_log_likelihood_function(data, tau+1, t+1)
-                tmp_cond = Pcp[j-1, tau] + P[tau+1, t] + Q[t + 1] + g[t - tau] - Q[tau + 1]
-                Pcp[j, t] = np.logaddexp(Pcp[j, t], tmp_cond)
+            tmp_cond = Pcp[j-1, j-1:t] + P[j:t+1, t] + Q[t + 1] + g[0:t-j+1] - Q[j:t+1]
+            Pcp[j, t] = logsumexp(tmp_cond)
 
     return Q, P, Pcp
 
 @dynamicProgramming
 def gaussian_obs_log_likelihood(data, t, s):
+    s += 1
     n = s - t
     mean = data[t:s].sum(0) / n
 
@@ -104,8 +104,8 @@ def gaussian_obs_log_likelihood(data, t, s):
     betaT = 1 + 0.5 * ((data[t:s] - mean) ** 2).sum(0) + ((n)/(1 + n)) * (mean**2 / 2)
     scale = (betaT*(nuT + 1))/(alphaT * nuT)
 
-    # splitting the PDF of the student distribution up is /much/ faster. (~ factor 20)
-    # using sum over for loop is even more worthwhile
+    # splitting the PDF of the student distribution up is /much/ faster.
+    # (~ factor 20) using sum over for loop is even more worthwhile
     prob = np.sum(np.log(1 + (data[t:s] - muT)**2/(nuT * scale)))
     lgA = gammaln((nuT + 1) / 2) - np.log(np.sqrt(np.pi * nuT * scale)) - gammaln(nuT/2)
 
@@ -115,8 +115,10 @@ def gaussian_obs_log_likelihood(data, t, s):
 def const_prior(r, l):
     return 1/(l)
 
+
 def geometric_prior(t, p):
     return p * ((1 - p) ** (t - 1))
 
+
 def neg_binominal_prior(t, k, p):
     return comb(t - k, k - 1) * p ** k * (1 - p) ** (t - k)

Example Code.ipynb

@@ -0,0 +1,14 @@
+#!/usr/bin/env python2
+
+from distutils.core import setup
+import bayesian_changepoint_detection
+
+setup(name='bayesian_changepoint_detection',
+      version=bayesian_changepoint_detection.__version__,
+      description='Some Bayesian changepoint detection algorithms',
+      author='Johannes Kulick',
+      author_email='johannes.kulick@ipvs.uni-stuttgart.de',
+      url='http://github.com/hildensia/bayesian_changepoint_detection',
+      packages = ['bayesian_changepoint_detection'],
+      requires=['scipy', 'numpy']
+     )