Using Matrix operation within pm_Model not working

Questions v5
1

For a changepoint time series analysis I was checking a modelling technique on PYMC3. Timeseers.

The Basic idea is to have an indicator matrix A of dimension: “Length of the time series” x “changepoints”.

I created an artificial time series with two changepoint after 1/3 and 2/3 of the time series
image

with pm.Model() as ts_model:
    #the initial growth and innteresection is k and m respectively  
    k = pm.Normal('k', mu=0, sigma=1)
    m = pm.Normal('m', mu=0, sigma=1)
    # after each change point  the growth rate  may change with delta 
    delta = pm.Laplace('delta', mu=0, b=0.01, shape=n_changepoints)
    # reflect the slope with new intersections
    gamma = delta * -s

    growth = k + pm.math.dot(A, delta)
    offset = m + pm.math.dot(A, gamma)
    
    trend = growth * t + offset

    error = pm.HalfCauchy('sigma', beta=0.5)
    pm.Normal('observation', mu=trend, sigma=error, observed=trend_rnd)
    idata = pm.sample(tune=1000)

The model runs fine but it doesn’t detect any changepoint. The delta array should react on the change in direction of the time series. But it doesn’t. I have the feeling that the matrix operation is somehow applied wrongly by me.

Any hint appreciated.

2

Can you provide a runnable example? It’s difficult to figure out what’s going on otherwise.

3

This code should show the issue.

import numpy as np
import pymc as pm
import matplotlib.pyplot as plt
import arviz as az

n_changepoints = 2

t = np.arange(50) / 50
s = np.linspace(0,max(t),n_changepoints+2)[1:-1]

#indicator matrix 
A = (t[:, None] > s) *1.0

k_det = 1
m_det = 0
delta_det = [-2.0, 1.]
gamma_det = delta_det * -s

g_det = (k_det + A@delta_det)
o_det = (m_det + A@gamma_det)


trend_det = g_det*t+o_det
trend_rnd = trend_det +  np.random.normal(0, 0.02, trend_det.size)
trend_rnd = (trend_rnd - np.mean(trend_rnd))/np.std(trend_rnd)
plt.plot(trend_rnd)

with pm.Model() as ts_model:
    #the initial growth and innteresection is k and m respectively  
    k = pm.Normal('k', mu=0, sigma=1)
    m = pm.Normal('m', mu=0, sigma=1)
    # after each change point  the growth rate  may change with delta 
    delta = pm.Laplace('delta', mu=0, b=0.01, shape=n_changepoints)
    # reflect the slope with new intersections
    gamma = delta * -s

    growth = k + pm.math.dot(A, delta)
    offset = m + pm.math.dot(A, gamma)
    
    trend = growth * t + offset

    error = pm.HalfCauchy('sigma', beta=0.5)
    pm.Normal('observation', mu=trend, sigma=error, observed=trend_rnd)
    idata = pm.sample(tune=1000)

az.plot_trace(idata,combined=True)
az.plot_forest(idata,combined=True)
4

I suspect that it’s because your priors are strongly informative. If you loosen up the b on the delta parameters, things seem a bit better. Not sure if that’s the only problem, but at least the posterior predictive samples look better at that point (and the posterior delta values are far from the prior`).

5

Hi Clumann, thanks a lot your feedback helped me already!

1 Like