quick_switch_point_finder.py

import pymc as pm
import numpy as np
import math

site_id = 911

def MCMC_switch_point_finder(count_data):
    n_count_data = len(count_data)
    alpha = 1.0 / count_data.mean()  # Recall count_data is the
                                   # variable that holds our txt counts
    lambda_1 = pm.Exponential("lambda_1", alpha)
    lambda_2 = pm.Exponential("lambda_2", alpha)
    tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data)
    @pm.deterministic
    def lambda_(tau=tau, lambda_1=lambda_1, lambda_2=lambda_2):
        out = np.zeros(n_count_data)
        out[:tau] = lambda_1  # lambda before tau is lambda1
        out[tau:] = lambda_2  # lambda after (and including) tau is lambda2
        return out
    observation = pm.Poisson("obs", lambda_, value=count_data, observed=True)
    model = pm.Model([observation, lambda_1, lambda_2, tau])
    mcmc = pm.MCMC(model)
    mcmc.sample(40000, 10000, 1)
    lambda_1_samples = mcmc.trace('lambda_1')[:]
    lambda_2_samples = mcmc.trace('lambda_2')[:]
    tau_samples = mcmc.trace('tau')[:]
    N = tau_samples.shape[0]
    expected_enters_per_day = np.zeros(n_count_data)
    for day in range(0, n_count_data):
        # ix is a bool index of all tau samples corresponding to
        # the switchpoint occurring prior to value of 'day'
        ix = day < tau_samples
        # Each posterior sample corresponds to a value for tau.
        # for each day, that value of tau indicates whether we're "before"
        # (in the lambda1 "regime") or
        #  "after" (in the lambda2 "regime") the switchpoint.
        # by taking the posterior sample of lambda1/2 accordingly, we can average
        # over all samples to get an expected value for lambda on that day.
        # As explained, the "enters count" random variable is Poisson distributed,
        # and therefore lambda (the poisson parameter) is the expected value of
        # "enters count".
        expected_enters_per_day[day] = (lambda_1_samples[ix].sum()
                                       + lambda_2_samples[~ix].sum()) / N
    return expected_enters_per_day, tau_samples


def quick_changepoint_finder(d):
    n = len(d)
    # dbar = sum(d)/float(n)
    dbar = np.mean(d)
    # dsbar = sum (d*d)/float(n)
    dsbar = np.mean(np.multiply(d,d))
    fac = dsbar-np.square(dbar)
    summ = 0
    summup = []
    for z in range(n):
        summ+=d[z]
        summup.append(summ)
    y = []
    for m in range(n-1):
        pos=m+1
        mscale = 4*(pos)*(n-pos)
        Q = summup[m]-(summ-summup[m])
        U = -np.square(dbar*(n-2*pos) + Q)/float(mscale) + fac
        y.append(-(n/float(2)-1)*math.log(n*U/2) - 0.5*math.log((pos*(n-pos))))
    z, zz = np.max(y), np.argmax(y)
    mean1 = sum(d[:zz+1])/float(len(d[:zz+1]))
    mean2=sum(d[(zz+1):n])/float(n-1-zz)
    return y, zz, mean1, mean2