Using Pymc3 to do forecasting and numerical integration

Best way I would say is to use pm.Data, which is design to handle prediction better:

training_set = df_res[~df_res.index.isin(t_pred_idx)]
test_set = df_res.loc[t_pred_idx]
with pm.Model() as model_logistic:
    predictor = pm.Data('predictor', training_set.t.values)
    observed = pm.Data('observed', training_set.y.values) 
    sigma =  pm.HalfNormal('sigma', sigma=5) # Half Normal are already bounded
    alpha = pm.Bound(pm.Normal, lower=0.0)('alpha', mu=1, sigma=1,shape=len(np.unique(indexer))) # pm.Normal('alpha',mu=1, sigma=1) # 1,1  # make boundary: https://docs.pymc.io/api/bounds.html
    beta =  pm.Bound(pm.Normal, lower=0.0)('beta', mu=1, sigma=1,shape=len(np.unique(indexer))) #pm.Normal('beta',mu=1, sigma=1) # 1,1
    t0 = pm.Bound(pm.Normal, lower=0.0)('t0', mu=10, sigma=10,shape=len(np.unique(indexer))) # pm.Normal('t0',mu=10, sigma=10) #normal(10,10);
    
    mu = alpha / (1 + pm.math.exp(-(beta*(predictor-t0))))
    # do indexing here instead
    Y_obs = pm.Normal('Y_obs', mu=mu[indexr], sigma=sigma, observed=observed)

The nice thing here is that, by pushing the indexing to observed, for prediction you can just use mu:

with model_logistic:
    pm.set_data({'predictor': test_set.t.values})
    Y_predict = pm.Normal('Y_predict', mu=mu, sigma=sigma, shape=...)

As for numerical integration, you can take a look at Custom theano Op to do numerical integration

3 Likes