Merge pull request hildensia#5 from hildensia/mixture_observations

Mixture observations

Author: hildensia

bayesian_changepoint_detection/cy_offline.pyx

@@ -0,0 +1,45 @@
+from __future__ import division
+import numpy as np
+
+def generate_normal_time_series(num, minl=50, maxl=1000):
+  data = np.array([], dtype=np.float64)
+  partition = np.random.randint(minl, maxl, num)
+  for p in partition:
+    mean = np.random.randn()*10
+    var = np.random.randn()*1
+    if var < 0:
+      var = var * -1
+    tdata = np.random.normal(mean, var, p)
+    data = np.concatenate((data, tdata))
+  return partition, np.atleast_2d(data).T
+
+def generate_multinormal_time_series(num, dim, minl=50, maxl=1000):
+  data = np.empty((1,dim), dtype=np.float64)
+  partition = np.random.randint(minl, maxl, num)
+  for p in partition:
+    mean = np.random.standard_normal(dim)*10
+    # Generate a random SPD matrix
+    A = np.random.standard_normal((dim,dim))
+    var = np.dot(A,A.T)
+
+    tdata = np.random.multivariate_normal(mean, var, p)
+    data = np.concatenate((data, tdata))
+  return partition, data[1:,:]
+
+def generate_xuan_motivating_example(minl=50, maxl=1000):
+  dim = 2
+  num = 3
+  partition = np.random.randint(minl, maxl, num)
+  mu = np.zeros(dim)
+  Sigma1 = np.asarray([[1.0,0.75],[0.75,1.0]])
+  data = np.random.multivariate_normal(mu, Sigma1, partition[0])
+  Sigma2 = np.asarray([[1.0,0.0],[0.0,1.0]])
+  data = np.concatenate((data,np.random.multivariate_normal(mu, Sigma2, partition[1])))
+  Sigma3 = np.asarray([[1.0,-0.75],[-0.75,1.0]])
+  data = np.concatenate((data,np.random.multivariate_normal(mu, Sigma3, partition[2])))
+  return partition, data
+
+
+
+
+

bayesian_changepoint_detection/generate_data.py

@@ -1,6 +1,6 @@
 from __future__ import division
 import numpy as np
-from scipy.special import gammaln
+from scipy.special import gammaln, multigammaln
 from scipy.misc import comb
 from decorator import decorator
 
@@ -39,12 +39,13 @@ def offline_changepoint_detection(data, prior_func,
     """Compute the likelihood of changepoints on data.
 
     Keyword arguments:
-    data -- the time series data
-    prior_func -- a function given the likelihood of a changepoint given the
-                  distance to the last one
+    data                                -- the time series data
+    prior_func                          -- a function given the likelihood of a changepoint given the distance to the last one
     observation_log_likelihood_function -- a function giving the log likelihood
                                            of a data part
-    P -- the likelihoods if pre-computed
+    truncate                            -- the cutoff probability 10^truncate to stop computation for that changepoint log likelihood
+
+    P                                   -- the likelihoods if pre-computed
     """
 
     n = len(data)
@@ -65,9 +66,9 @@ def offline_changepoint_detection(data, prior_func,
     Q[n-1] = P[n-1, n-1]
 
     for t in reversed(range(n-1)):
-        P_next_cp = -np.inf  # == -log(0)
+        P_next_cp = -np.inf  # == log(0)
         for s in range(t, n-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]
@@ -82,7 +83,7 @@ def offline_changepoint_detection(data, prior_func,
 
         # (1 - G) is numerical stable until G becomes numerically 1
         if G[n-1-t] < -1e-15:  # exp(-1e-15) = .99999...
-            antiG = np.log(1 - np.exp(G[n-1-t])) 
+            antiG = np.log(1 - np.exp(G[n-1-t]))
         else:
             # (1 - G) is approx. -log(G) for G close to 1
             antiG = np.log(-G[n-1-t])
@@ -108,28 +109,65 @@ def gaussian_obs_log_likelihood(data, t, s):
     s += 1
     n = s - t
     mean = data[t:s].sum(0) / n
-    
+
     muT = (n * mean) / (1 + n)
     nuT = 1 + n
     alphaT = 1 + n / 2
     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
     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)
-    
+
     return np.sum(n * lgA - (nuT + 1)/2 * prob)
 
+def ifm_obs_log_likelihood(data, t, s):
+    '''Independent Features model from xuan et al'''
+    s += 1
+    n = s - t
+    x = data[t:s]
+    if len(x.shape)==2:
+        d = x.shape[1]
+    else:
+        d = 1
+        x = np.atleast_2d(x).T
+
+    N0 = d          # weakest prior we can use to retain proper prior
+    V0 = np.var(x)
+    Vn = V0 + (x**2).sum(0)
+
+    # sum over dimension and return (section 3.1 from Xuan paper):
+    return d*( -(n/2)*np.log(np.pi) + (N0/2)*np.log(V0) - \
+        gammaln(N0/2) + gammaln((N0+n)/2) ) - \
+        ( ((N0+n)/2)*np.log(Vn) ).sum(0)
+
+def fullcov_obs_log_likelihood(data, t, s):
+    '''Full Covariance model from xuan et al'''
+    s += 1
+    n = s - t
+    x = data[t:s]
+    if len(x.shape)==2:
+        dim = x.shape[1]
+    else:
+        dim = 1
+        x = np.atleast_2d(x).T
+
+    N0 = dim          # weakest prior we can use to retain proper prior
+    V0 = np.var(x)*np.eye(dim)
+    Vn = V0 + np.array([np.outer(x[i], x[i].T) for i in xrange(x.shape[0])]).sum(0)
+
+    # section 3.2 from Xuan paper:
+    return -(dim*n/2)*np.log(np.pi) + (N0/2)*np.linalg.slogdet(V0)[1] - \
+        multigammaln(N0/2,dim) + multigammaln((N0+n)/2,dim) - \
+        ((N0+n)/2)*np.linalg.slogdet(Vn)[1]
 
 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)

bayesian_changepoint_detection/offline_changepoint_detection.py

@@ -0,0 +1,43 @@
+from __future__ import division
+import numpy as np
+import matplotlib.pyplot as plt
+import seaborn
+
+import cProfile
+import bayesian_changepoint_detection.offline_changepoint_detection as offcd
+import bayesian_changepoint_detection.generate_data as gd
+from functools import partial
+
+if __name__ == '__main__':
+  show_plot = True
+  dim = 4
+  if dim == 1:
+    partition, data = gd.generate_normal_time_series(7, 50, 200)
+  else:
+    partition, data = gd.generate_multinormal_time_series(7, dim, 50, 200)
+  changes = np.cumsum(partition)
+
+  if show_plot:
+    fig, ax = plt.subplots(figsize=[16,12])
+    for p in changes:
+      ax.plot([p,p],[np.min(data),np.max(data)],'r')
+    for d in range(dim):
+      ax.plot(data[:,d])
+    plt.show()
+
+
+  #Q, P, Pcp = offcd.offline_changepoint_detection(data,partial(offcd.const_prior, l=(len(data)+1)),offcd.gaussian_obs_log_likelihood, truncate=-20)
+  #Q_ifm, P_ifm, Pcp_ifm = offcd.offline_changepoint_detection(data,partial(offcd.const_prior, l=(len(data)+1)),offcd.ifm_obs_log_likelihood,truncate=-20)
+  Q_full, P_full, Pcp_full = offcd.offline_changepoint_detection(data,partial(offcd.const_prior, l=(len(data)+1)),offcd.fullcov_obs_log_likelihood, truncate=-50)
+
+  if show_plot:
+    fig, ax = plt.subplots(figsize=[18, 16])
+    ax = fig.add_subplot(2, 1, 1)
+    for p in changes:
+      ax.plot([p,p],[np.min(data),np.max(data)],'r')
+    for d in range(dim):
+      ax.plot(data[:,d])
+    ax = fig.add_subplot(2, 1, 2, sharex=ax)
+    ax.plot(np.exp(Pcp_full).sum(0))
+    plt.show()
+

example.py

@@ -0,0 +1,34 @@
+''' Example from Xiang Xuan's thesis: Section 3.2'''
+from __future__ import division
+import numpy as np
+import matplotlib.pyplot as plt
+
+import bayesian_changepoint_detection.offline_changepoint_detection as offcd
+import bayesian_changepoint_detection.generate_data as gd
+from functools import partial
+
+if __name__ == '__main__':
+  show_plot = True
+
+  partition, data = gd.generate_xuan_motivating_example(50,200)
+  changes = np.cumsum(partition)
+
+  Q_ifm, P_ifm, Pcp_ifm = offcd.offline_changepoint_detection(data,partial(offcd.const_prior, l=(len(data)+1)),offcd.ifm_obs_log_likelihood,truncate=-20)
+  Q_full, P_full, Pcp_full = offcd.offline_changepoint_detection(data,partial(offcd.const_prior, l=(len(data)+1)),offcd.fullcov_obs_log_likelihood, truncate=-20)
+
+  if show_plot:
+    fig, ax = plt.subplots(figsize=[18, 16])
+    ax = fig.add_subplot(3, 1, 1)
+    for p in changes:
+      ax.plot([p,p],[np.min(data),np.max(data)],'r')
+    for d in range(2):
+      ax.plot(data[:,d])
+    plt.legend(['Raw data with Original Changepoints'])
+    ax1 = fig.add_subplot(3, 1, 2, sharex=ax)
+    ax1.plot(np.exp(Pcp_ifm).sum(0))
+    plt.legend(['Independent Factor Model'])
+    ax2 = fig.add_subplot(3, 1, 3, sharex=ax)
+    ax2.plot(np.exp(Pcp_full).sum(0))
+    plt.legend(['Full Covariance Model'])
+    plt.show()
+