Switch point Metropolis tuning


#1

I have data with potentially 3 switch points:
data
The model will automatically choose Metropolis for the sampling of the tau_1 and tau_2, and NUTS for the lambdas.

N = np.arange(0, all_points)
with pm.Model() as model:
alpha = 1.0/data.mean()
lambda_1 = pm.Exponential(“lambda_1”, alpha)
lambda_2 = pm.Exponential(“lambda_2”, alpha)
lambda_3 = pm.Exponential(“lambda_3”, alpha)

tau1 = pm.DiscreteUniform("tau1", lower=N.min(), upper=N.max()) 
tau2 = pm.DiscreteUniform("tau2", lower=tau1, upper=N.max()) 

_mu = T.switch(tau1>=N,lambda_1,lambda_2)
mu = T.switch(tau2>=N,_mu,lambda_3)

observation = pm.Poisson("obs", mu, observed=data)
trace = pm.sample(10000, tune=1000)

When tune=500 it doesn’t converge, but with tune=1000, both tau’s converge successfully.

Can you please share with me and the Pymc3 community how tuning works or what should be considered when tuning a Metropolis. Do you know of articles/websites where I could read more about tuning of the Metropolis?


#2

I am surprised that discrete switch point works well here, still, I suggest you to use a continuous step function instead of switch, more details see: https://stackoverflow.com/questions/49144144/convert-numpy-function-to-theano/49152694#49152694

Back to your original question, in pm.sample each sampler has some target that it would try to reach during tuning. For example, in Metropolis, it adjusts the scale of the transition kernel so that the acceptance probability is around 50%

In NUTS, tuning is much more complex as it optimizes the step size for leapfrog integrator, the mass matrix of the Hamiltonian transition kernel. The best reference for this is the Stan User Manuel http://mc-stan.org/users/documentation/ chapter 34.2, although we dont use the excat tuning algorithm as Stan.


#3

When I try to run what you created:

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

with pm.Model() as m2:
    lambda0 = pm.Normal('lambda0', mu, sd=sd)
    lambdad = pm.Normal('lambdad', 0, sd=sd, shape=nbreak-1)
    trafo = Composed(pm.distributions.transforms.LogOdds(), Ordered())
    b = pm.Beta('b', 1., 1., shape=nbreak-1, transform=trafo,
                    testval=[0.3, 0.5])
    theta_ = pm.Deterministic('theta', tt.exp(lambda0 +
                                              logistic(lambdad[0], b[0]) +
                                              logistic(lambdad[1], b[1])))
    obs = pm.Poisson('obs', theta_, observed=y)
    
    trace = pm.sample(1000, tune=1000)

I get the following error:

AttributeError: Can’t get attribute ‘Composed’ on <module ‘main’ (built-in)>

Is there any way around that?


#4

Yes, more details is in the notebook, link copied below:


#5

I’m sorry for the confusion. I already ran the whole notebook the first time, but it gives the error as soon as you try to run the model.

AttributeError: Can’t get attribute ‘Composed’ on <module ‘main’ (built-in)>

I pretty much copy-pasted everything but the error keeps repeating. I tried changing some parameters but I feel the problem is within the class Composed(Transform):

Not sure…


#6

That’s pretty strange… what is your pymc3 version?
Also, could you please paste the full error trace?


#7

Pymc3 version 3.3
I’m running it in jupyter notebook (version 5.3.1).
When I run the model, it gets up to this part:

Auto-assigning NUTS sampler…
Initializing NUTS using jitter+adapt_diag…
c:\anaconda\lib\site-packages\pymc3\model.py:384: FutureWarning: Conversion of the second argument of issubdtype from float to np.floating is deprecated. In future, it will be treated as np.float64 == np.dtype(float).type.
if not np.issubdtype(var.dtype, float):
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [b_logodds_ordered__, lambdas]

Then it just gets stuck, but if I you look at the Jupyter command prompt you’d see:

Process SpawnPoolWorker-2:
Traceback (most recent call last):
File “c:\anaconda\lib\multiprocessing\process.py”, line 249, in _bootstrap
self.run()
File “c:\anaconda\lib\multiprocessing\process.py”, line 93, in run
self._target(*self._args, **self._kwargs)
File “c:\anaconda\lib\multiprocessing\pool.py”, line 108, in worker
task = get()
File “c:\anaconda\lib\site-packages\joblib\pool.py”, line 362, in get
return recv()
File “c:\anaconda\lib\multiprocessing\connection.py”, line 251, in recv
return _ForkingPickler.loads(buf.getbuffer())
AttributeError: Can’t get attribute ‘Composed’ on <module ‘main’ (built-in)>

Do you recommend I update Jupyter or Pymc3 versions?


#8

Looks like you have problem of the multiprocessing - try upgrading joblib.
For now, setting njobs=1 in pm.sample() should temporally suppress the error.


#9

Thank you, junpenglao, and thank you for your patience. Upgrading joblib didn’t do it, but it worked when changed the njobs=1.
When I try to apply your model to other datasets that I am analyzing, is it correct that the testval would most likely change for ‘b’?

In the current model, the testval for ‘b’ is [0.3,0.5]:

with pm.Model() as m2:
    lambda0 = pm.Normal('lambda0', mu, sd=sd)
    lambdad = pm.Normal('lambdad', 0, sd=sd, shape=nbreak-1)
    trafo = Composed(pm.distributions.transforms.LogOdds(), Ordered())
    b = pm.Beta('b', 1., 1., shape=nbreak-1, transform=trafo,
                    testval=[0.3, 0.5])

Can you please share why you chose testval=[0.3, 0.5]. I know these are starting values for ‘b’, but would those change if my initial observed data was different?


#10

The testval here is to avoid a bug in the Ordered transformed: the default testval of the Beta(1, 1) is .5, and with shape=2 it would be [.5, .5], however it generate an error for Ordered as the two value is the same. So as long as you specify a value that is sorted but not the same then it would work.


#11

thank you. :slight_smile: