Switch point Metropolis tuning

Questions
12

Jutting in here as I have a very similar model that I am working with. I was previously using the DiscreteUniform distribution to model the switchpoints in my data, but now am trying to shift to the logistic as explained here.

I am trying to model my data with a single regression discontinuity - so two straight lines fit my data, and at some point in time, the model has to switch from one straight line to the other. Here’s an example:
image

I tried the same logistic idea as in @junpenglao 's gist, here’s the code:

    def logistic(L, x0, k=500, t=np.linspace(0., 1., 2500)):
         return L/(1+tt.exp(-k*(t-x0)))

And the model:

    In [34]: with pm.Model() as model_before_700:                    
        ...:     alpha = pm.Normal("alpha", mu = 0, sd = 3.0, shape = 2)           
        ...:     beta = pm.Normal("beta", mu = 0, sd = 1.0, shape = 2)
        ...:     switchpoints = pm.Beta("switchpoints", 1, 1)                       
        ...:     sd = pm.HalfCauchy("sd", 0.5, shape = 2)
        ...:     intercept = alpha[0] + logistic(alpha[1], switchpoints)
        ...:     slope = beta[0] + logistic(beta[1], switchpoints)
        ...:     dev = sd[0] + logistic(sd[1], switchpoints)
        ...:     regression = intercept + slope*np.arange(2500)                    
        ...:     observed = pm.Normal("observed", mu = regression, sd = dev, observed = data_before_700)                
    In [34]: with model_before_700:                                  
        ...:     trace_before_700 = pm.sample(tune = 2000, draws = 1000, njobs = 3)
        ...:

I consistently get:
ValueError: Mass matrix contains zeros on the diagonal. Some derivatives might always be zero.

I have seen previous discussions saying that this might be due to overflow errors, so I think I might be missing something simple in my model. Any help is much appreciated :slight_smile:

Here’s relevant versions etc:

    In [38]: pm.__version__
    Out[38]: '3.3'

    In [39]: theano.__version__
    Out[39]: '1.0.1'

    In [40]: joblib.__version__
    Out[40]: '0.11'