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