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?