# Help with Change Point Analysis

Hi everyone, new to Julia and Turing.jl here! I have a question about why my change point model will not converge properly. This question was asked before here, but none of the answers there worked either.

I am creating some synthetic data using 2 Poisson distributions and merging them together into 1 long list with a change point halfway through. See line 4 of code for how the data are defined / created (I visually inspected the scatter plot and the data appear to be generated correctly). I am then defining a model that has a prior for lambda_1, lambda_2, and tau (the change point).

``````using Random, Distributions
Random.seed!(24)
# creating 200 synthetic data points with a change point halfway through
data = vcat(rand(Poisson(20), 100), rand(Poisson(2), 100))

# defining the model
using Turing
@model function cp_model(y)
λ₁ ~ Gamma(0.01, 100)
λ₂ ~ Gamma(0.01, 100)
τ  ~ DiscreteUniform(0,200)

for idx in length(y)
λ = idx < τ ? λ₁ : λ₂
y[idx] ~ Poisson(λ)
end
return λ₁,λ₂,τ
end

# sampling using metropolis hastings
chains = sample(cp_model(data), MH(), init_params = (λ₁ = 12, λ₂ = 2, τ = 100), 10000)
``````

Everything is running just fine, but the results are SUPER wonky. Both lambdas have a mean less than 1 (small standard deviations as well), but tau is somehow right on the money at 101. Can anyone help me out with this?

Welcome @GrudgingDesert1! First, I notice your model contains a bug: you want `1:length(y)` in the loop (or even better, `eachindex(y)`). Even after fixing that, though, I still get poor convergence.

In general, I’ve found the `MH()` sampler in Turing challenging to use, since it doesn’t do automatic step-size tuning. The fix I’d suggest here is to just make `τ` a continuous variable and sample using `NUTS()`:

``````@model function cp_model2(y)
λ₁ ~ Gamma(0.01, 100)
λ₂ ~ Gamma(0.01, 100)
τ  ~ Uniform(0,200)

for idx in eachindex(y)
λ = idx < τ ? λ₁ : λ₂
y[idx] ~ Poisson(λ)
end
return λ₁,λ₂,τ
end

chains2 = sample(cp_model2(data), NUTS(), 1000)
``````

For me this converges well and recovers all the parameters.

1 Like

Thank you for that! I didn’t realize how bad the metropolis was performing! replacing tau with a continuous uniform prior and using NUTS did the trick. I did find another way to make this work that leaves the tau prior as a discrete variable. I used a Gibbs sampler: NUTS sampling the lambdas and metropolis sampling tau. Here is the code for that:

``````using Random, Distributions
Random.seed!(24)
# creating 200 synthetic data points with a change point halfway through
data = vcat(rand(Poisson(20), 100), rand(Poisson(2), 100))

# defining the model
using Turing
@model function model(y)
λ₁ ~ Gamma(0.001, 1000)
λ₂ ~ Gamma(0.001, 1000)
τ  ~ DiscreteUniform(1, 200)

for idx in eachindex(y)
λ = idx < τ ? λ₁ : λ₂
y[idx] ~ Poisson(λ)
end
return λ₁,λ₂,τ
end

# sampling using NUTS
chains = sample(model(data),
Gibbs(NUTS(1000, 0.65, :λ₁, :λ₂), MH(:τ)),
init_params = (λ₁ = 20, λ₂ = 2, τ = 100),
10000)
``````

Both of them yielding essentially the same results with the added benefit of leaving tau as a discrete variable. Thanks of you help on this!

1 Like