pymc

change point detection with two transitions

this question was asked previously and not answered (5 June) but maybe putting it in context makes more sense. I have done the change point tutorial with the two lambdas and extended with 2 change point so the modelling is now:

# the exp parameter expected is the inverse of the average from sampled series
alpha = 1.0 / count_data.mean() 
# regime 1 poisson
lambda_1 = pm.Exponential("lambda_1", alpha)
# regime 2 poisson
lambda_2 = pm.Exponential("lambda_2", alpha)

# regime 3 poisson
lambda_3 = pm.Exponential("lambda_3", alpha)

# change point is somewhere in between with equal probabilities
tau1 = pm.DiscreteUniform("tau1", lower=0, upper=n_count_data)
# change point is somewhere in between with equal probabilities
tau2 = pm.DiscreteUniform("tau2", lower=0, upper=n_count_data)

@pm.deterministic
def lambda_(tau1=tau1,tau2=tau2, lambda_1=lambda_1, lambda_2=lambda_2):
    out = np.zeros(n_count_data)
    out[:tau1] = lambda_1  # lambda before tau is lambda1
    out[tau1:tau2] = lambda_2  # lambda between periods is lambda2
    out[tau2:] = lambda_3  # lambda after (and including) tau2 is lambda3
    return out

observation = pm.Poisson("obs", lambda_, value=count_data, observed=True)

model = pm.Model([observation, lambda_1, lambda_2, tau1,tau2])

# markov monte carlo chain
mcmc = pm.MCMC(model)
mcmc.sample(40000, 10000, 1)

The question is that in the deterministic variable how do I actually tell the model that I only need to consider when tau1 is less than tau2? The problem is that when tau2 precedes tau1 there is a time symmetry which is computationally non necessary.

Any help is welcome.

Solution

  • I haven't tested it, but I think you could do something like this:

    # change point is somewhere in between with equal probabilities
tau1 = pm.DiscreteUniform("tau1", lower=0, upper=n_count_data)
# change point is somewhere in between with equal probabilities
tau2 = pm.DiscreteUniform("tau2", lower=tau1, upper=n_count_data)

    That way tau2 is constrained to be at least as large as tau1. You may have to think a little bit about whether tau1 and tau2 should be allowed to coincide.