# Help with Change Point Analysis

**URL:** <https://discourse.julialang.org/t/help-with-change-point-analysis/113677>\
**Category:** New to Julia\
**Tags:** question, turing\
**Created:** [May 1, 2024, 12:06am UTC](https://discourse.julialang.org/t/help-with-change-point-analysis/113677 "2024-05-01T00:06:16Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![GrudgingDesert1](https://avatars.discourse-cdn.com/v4/letter/g/aca169/32.png) [@GrudgingDesert1](https://discourse.julialang.org/u/GrudgingDesert1)\
**Post date:** [May 1, 2024, 12:06am UTC](https://discourse.julialang.org/t/help-with-change-point-analysis/113677/1 "2024-05-01T00:06:16Z")

</div>

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](https://discourse.julialang.org/t/help-converting-pymc3-switchpoint-analysis-code-to-turing-jl/39880), 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).

```julia
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?

---

<div class="post-metadata">

**Author:** ![ElOceanografo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eloceanografo/32/624_2.png) [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)\
**Post date:** [May 1, 2024, 12:36am UTC](https://discourse.julialang.org/t/help-with-change-point-analysis/113677/2 "2024-05-01T00:36:37Z")

</div>

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()`:

```julia
@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.

---

<div class="post-metadata">

**Author:** ![GrudgingDesert1](https://avatars.discourse-cdn.com/v4/letter/g/aca169/32.png) [@GrudgingDesert1](https://discourse.julialang.org/u/GrudgingDesert1)\
**Post date:** [May 1, 2024, 4:36am UTC](https://discourse.julialang.org/t/help-with-change-point-analysis/113677/3 "2024-05-01T04:36:52Z")

</div>

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:

```julia
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!