change point detection with two transitions

Question

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.

score 1 · Accepted Answer · answered Aug 03 '14 at 20:01

1

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.

answered Aug 03 '14 at 20:01

jcrudy

3,921
1
24
31

Excellent that worked quite well under a simplified assumption where I have set the gap to be a constant of 1 day but I will also try to model that as a uniform or gaussian random variable. Full model is on my own answer. – user2471214 Aug 08 '14 at 14:00

score 0 · Answer 2 · answered Aug 08 '14 at 14:05

The full model under the assumption of a deterministic gap between the taus follows:

# 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=tau1+1, upper=n_count_data)

@pm.deterministic
def lambda_(tau1=tau1,tau2=tau2, lambda_1=lambda_1, lambda_2=lambda_2,lambda_3=lambda_3):
    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,lambda_3, tau1,tau2])

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

lambda_1_samples = mcmc.trace('lambda_1')[:]
lambda_2_samples = mcmc.trace('lambda_2')[:]
lambda_3_samples = mcmc.trace('lambda_3')[:]

tau1_samples = mcmc.trace('tau1')[:]
tau2_samples = mcmc.trace('tau2')[:]

Will also try with the random gap and see how it goes.

Jonas Lindeløv · Answer 3 · 2020-01-10T10:36:41.713

If you are open to use R to solve the same inference problem, the mcp package provides a higher-level interface for change point problems. It has order-restricted change point parameters by default.

Here is a model for three intercepts (two change points)

model = list(
  count ~ 1,
  ~ 1,
  ~ 1
)

library(mcp)
fit = mcp(model, data, family = poisson())

More info:

about Poisson models in mcp.
priors in mcp contains more about finer control of order-restriction.

change point detection with two transitions

3 Answers3