Help with Change Point Analysis

New to Julia
,
1

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?

2

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
3

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