具有两个转换的变化点检测

Question

这个问题之前被问过，但没有回答（6 月 5 日），但也许把它放在上下文中更有意义。 我已经用两个 lambda 完成了变化点教程，并用 2 个变化点进行了扩展，所以现在的建模是：

# 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)

问题是，在确定性变量中，我如何真正告诉模型我只需要在 tau1 小于 tau2 时考虑？ 问题是当 tau2 在 tau1 之前存在时间对称性，这在计算上是不必要的。

欢迎任何帮助。

Answer 1

我还没有测试过，但我认为你可以这样做：

# 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)

这样 tau2 被限制为至少与 tau1 一样大。您可能需要考虑一下是否应该允许 tau1 和 tau2 重合。

Answer 2

假设 taus 之间存在确定性差距的完整模型如下：

# 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')[:]

也会尝试随机间隙，看看效果如何。

Answer 3

如果您愿意使用 R 解决相同的推理问题，mcp 包为更改点问题提供了更高级别的接口。默认有顺序限制的变化点参数。

这里是三个截距（两个变化点）的模型

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

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

更多信息：

• about Poisson models in mcp。
• priors in mcp 包含更多关于更好地控制订单限制的内容。