Change Detection¶

from collections import defaultdict
np.random.seed(seed=111111)
np.set_printoptions(precision=3, suppress=True)

def dict_to_arrays(ddict):
    """
    Convenience function used by online simulator to bundle 
    residuals into a dict before returning
    """
    new_dict = {}
    for k,v in ddict.iteritems():
        new_dict[k] = np.array(v)
    return new_dict

def online_simulator(signal, change_detector): 
    """
    Function that simulates an online streaming scenario for change detection experiments.
    --- 
    Given a signal and a change detector, this simulator passes one signal data point at a time to
    the change detector and processes the results. 
    """
    #Initiate
    #change_detector = cd_algorithm()
    
    all_residuals = defaultdict(list)
    
    xx = 0
    #Iterate through the signal, passing data points to the algorithm. 
    for value in signal: 
        
        #calculate residuals, compare residuals with the stopping rule(s)
        check_results = next(change_detector.step(value))
        
        #process results
        rule_triggered    = check_results[0]
        res               = check_results[1]
        
        #store residuals
        for k,v in res.iteritems():
            all_residuals[k].append(v)
              
        if rule_triggered == True: 
            #stopping rule was triggered          
            return (True, dict_to_arrays(all_residuals))    

    #Rule wasn't triggered by end of signal
    return (False, dict_to_arrays(all_residuals))

def run_online_simulation(signal, change_detector, scale=True): 
    """Run simulation and print results"""
    
    #Run simulation
    results = online_simulator(signal, change_detector)
    
    #Display results
    print_sim_results(signal, results, scale=scale)
    
    #Return residuals
    residuals = results[1]
    return residuals

class change_detector(object):
    """
    A change detection algorithm. 
    
    The algorithm calculates residuals and updates them for each new value passed. 
    Residuals are checked against stopping rules at each change, yielding either True or False, accordingly. 
    
    """
    
    def __init__(self): 
        #Interim and calculated values
        self.signal_size = 0
        self.total_val = 0
        
        self.mean_ = np.nan
    
    def update_residuals(self, new_signal_value): 
        #Update residuals
        self.signal_size += 1
        self.total_val += new_signal_value
        self.mean_ = self.total_val / self.signal_size
    
    def _get_residual_dict(self): 
        """create a dictionary of residuals to return. 
        Inclues all class and instance variables ending in '_'
        """
        residuals_dict = {}
        for k,v in self.__dict__.iteritems():
            if k.endswith('_'):
                residuals_dict[k] = v
        
        return residuals_dict
        
    def check_stopping_rules(self, new_signal_value): 
        rules_triggered = False
        
        #No rules implemented
        pass
        
        return rules_triggered
    
    def _step(self, new_signal_value): 
        
        #update residuals
        self.update_residuals(new_signal_value)
        
        ## compare residuals to stopping_rules
        rules_triggered = self.check_stopping_rules(new_signal_value)
        
        if rules_triggered: 
            yield (True, self._get_residual_dict())
        
        else: 
            yield (False, self._get_residual_dict())
      
    def step(self, new_signal_value):
        return self._step(new_signal_value)

def print_sim_results(signal, results, **kwargs):
    """
    Another convenience function to print out the results of our experiment. 
    """
    #Get results
    stopped = results[0]
    residuals = results[1]
    print "Residuals: {}".format([res for res in residuals.viewkeys()])

    if stopped: 
        #length of residuals array tells us when the rule was triggered
        stop_point = len(residuals.itervalues().next())
          
        print "Change detected. Stopping Rule triggered at {}.\n".format(stop_point)
        plot_signal_and_residuals(signal, residuals, stop_point, **kwargs)
    else: 
        print "Stopping rule not triggered."
        plot_signal_and_residuals(signal, residuals, **kwargs)

def plot_signal_and_residuals(signal, residuals=None, stop_point=None, scale=True):
    """Convenience function to generate plots of the signal and the residuals"""
    
    if residuals is None:
        plotcount = 1
    else: 
        plotcount = 1 + len(residuals)
    
    fig, axes = plt.subplots(nrows=plotcount, 
                             ncols = 1, 
                             sharex=True,
                             figsize=(6, plotcount*3)
                             )

    #First plot the signal
    if plotcount > 1: 
        ax = axes[0]
    elif plotcount == 1: 
        ax = axes
        
    ax.plot(signal)
    ax.set_title('Signal')
    
    #Scale signal
    ax.set_ylim(signal.min()*.5, signal.max()*1.5)
    ax.set_xlim(0, len(signal))
        
    #Plot a horizontal line where the stop_point is indicated
    if stop_point is not None: 
        assert (stop_point > 0) & (stop_point < len(signal))
        ax.vlines(x=stop_point, ymin=0, ymax=ax.get_ylim()[1], 
                  colors='r', linestyles='dotted')
    
    #Now plot each residual
    if residuals is not None: 
        for ii, (res_name, res_values) in enumerate(residuals.iteritems()):
            ax = axes[ii+1]
            ax.plot(res_values)
            ax.set_title("Residual #{}: {}".format(ii+1, res_name))
            if scale: 
                ax.set_ylim(res_values.min()*0.5, res_values.max() * 1.5)
            ax.vlines(x=stop_point, ymin=0, ymax=ax.get_ylim()[1], 
                      colors='r', linestyles='dotted')
        

class cd_static_mean_detector(change_detector):
    """
    A change detection algorithm. 
    
    The algorithm calculates residuals and updates them for each new value passed. 
    Residuals are checked against stopping rules at each change, yielding either True or False, accordingly. 
    
    """
    
    def __init__(self, threshold=0.05): 
        #hyper-parameter(s)
        self.threshold = threshold
            
        #Interim and calculated values
        self.signal_size = 0
        self.total_val = 0
        
        #... and residuals
        self.diff_ = np.nan #np.zeros(1)
        self.mean_ = np.nan #np.zeros(1)
    
    def update_residuals(self, new_signal_value): 
        #Update residuals 
        self.signal_size += 1
        self.total_val += new_signal_value
        self.mean_ = self.total_val / self.signal_size
        self.diff_ = np.absolute(self.mean_ - new_signal_value)
    
    def check_stopping_rules(self, new_signal_value): 
        rules_triggered = False
        #check if new value is more than % different from mean
        threshold_level = self.mean_ * self.threshold
        
        if self.diff_ > threshold_level:
            #import pdb; pdb.set_trace();
            rules_triggered = True
        return rules_triggered
          
    def step(self, new_signal_value):
        return self._step(new_signal_value)

def welford(vals):
    #via https://gist.github.com/alexalemi/2151722
    
    m = vals[0]
    s = 0
    for k in range(1, len(vals)):
        oldm = m
        x = vals[k]
        newm = oldm + (x - oldm) / (k + 1)
        s = s + (x - newm) * (x - oldm)
        m = newm
    return s / (k+1)

class welford_demo(change_detector):
    """
    Testing Welford's algorithm
    implemented like our other change detector
    """
    
    def __init__(self): 
            
        #Interim and calculated values
        self.k = 0    #k is the signal_size

        #For Welford's 
        self.mean_ = 0.0     
        self.s = 0.0      # var = s / (k + 1)
                
    def update_residuals(self, new_signal_value): 
        #Update residuals 
        
        #Welford's. 
        oldm = self.mean_
        x = new_signal_value
        newm = oldm + (x - oldm) / (self.k + 1)
        self.s = self.s + (x - newm) * (x - oldm)
        self.mean_ = newm
        
        self.var_ = self.s / (self.k+1)
        
        self.k += 1
        

from collections import deque

class z_score_detector(change_detector):
    """
    A change detection algorithm based on z-scores
    ---
    1. We will use Welford's algorithm to collect mean and variance 
    of the full stream. 
    
    2. At the same time, we will collect a window of recent
    observations and compute local statistics on this sample. 
    
    3. Then, we compute a Z-score to compare the sample (window) 
    distribution to the population (global) distribution. 

    Our stopping rule is based on this computed z-score
    """
    
    def __init__(self, window_size = 100, threshold=0.05): 
        #hyper-parameter(s)
        self.threshold = threshold
        self.window_size = window_size
            
        #Interim and calculated values
        self.k = 0    #k is the signal_size
        
        #For Welford's 
        self.g_mean = 0.0     
        self.s = 0.0      # var = s / (k + 1)
        
        #For the Window
        self.window = deque(maxlen = window_size)
        
        #... and residuals
        self.z_score_ = np.nan
        
    def update_residuals(self, new_signal_value): 
        """update residuals """
        x = new_signal_value
        
        #Welford's on global window
        oldm = self.g_mean
        newm = oldm + (x - oldm) / (self.k + 1)
        self.s = self.s + (x - newm) * (x - oldm)
        self.g_mean = newm
        
        #Local window (sample)
        self.window.append(x)
        
        #Calculate z-score
        w_mean = np.mean(self.window)
        w_std = np.std(self.window)
        
        mean_diff = (self.g_mean - w_mean) / self.g_mean
        
        g_std = np.sqrt(self.s / (self.k+1))  #global std
        std_diff = (g_std - w_std) / g_std        
        SE = g_std / np.sqrt(self.window_size)

        self.z_score_ = (w_mean - self.g_mean) / SE
        self.window_mean_ = w_mean
        self.global_var_ = self.s / (self.k+1)
        
        self.k += 1
    
    def check_stopping_rules(self, new_signal_value): 
        rules_triggered = False

        #Stopping rule! 
        if np.abs(self.z_score_) > self.threshold:
            #import pdb; pdb.set_trace();
            rules_triggered = True
            
        return rules_triggered