Disaggregate and find the change point in time series using PyMC

Questions v3
2

Everything you’re doing is basically correct, you’re just not doing the changepoint correctly. You want to compare the value of tau to an array of dates, not to the correct date. Here’s a corrected model:

with pm.Model() as model:
    # Priors for change point
    t = pm.MutableData('t', np.arange(n-1))
    tau = pm.DiscreteUniform('tau', lower=0, upper=n - 1)
    
    # Priors for the means and standard deviations
    mu = pm.Normal("mean1", mu=4, sigma=2, size=2)
    sigma = pm.HalfNormal("var1", sigma=0.5, size=2)

    mu_obs = pm.math.switch(t < tau, mu[0], mu[1])
    sigma_obs = pm.math.switch(t < tau, sigma[0], sigma[1])
    obs = pm.Normal("obs", mu=mu_obs, sigma=sigma_obs, observed=true_data)

    # Sample from the posterior
    idata = pm.sample()

This gives the following results:

            mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_bulk  ess_tail  r_hat
tau       33.000  0.000  33.000   33.000      0.000    0.000    4000.0    4000.0    NaN
mean1[0]   6.313  0.156   6.021    6.610      0.002    0.001    5555.0    3010.0    1.0
mean1[1]   0.927  0.082   0.781    1.094      0.001    0.001    5172.0    2990.0    1.0
var1[0]    0.888  0.106   0.704    1.086      0.001    0.001    6786.0    3106.0    1.0
var1[1]    0.676  0.058   0.569    0.786      0.001    0.001    5354.0    3570.0    1.0

You can see that all the true parameters are recovered, and the change point is nailed down with zero uncertainty.