I have built a nice model with BART. I am trying to plot the “forecast” variable over time for the 4 forecast steps. I wanted to ask if the way I do plot the “forecast” is valid?
# BART :-D
forecast_steps = 4
with pm.Model() as V_Research:
σ = pm.HalfNormal("σ", Y.std())
μ = pmb.BART("μ", X_lagged, Y, m=50)
y = pm.Normal("y", μ, σ, observed=Y)
forecast = pm.Normal("forecast", mu=μ[-1], shape=(forecast_steps,))
idata_V_Research = pm.sample( return_inferencedata=True)
# Get the mean forecast
mean_forecast = np.zeros(forecast_steps)
for i in range(forecast_steps):
mean_forecast[i] = np.mean(idata_V_Research.posterior["forecast"][:, i])
# Plot the mean forecast
plt.plot(mean_forecast)
plt.xlabel("Forecast step")
plt.ylabel("Mean forecast")
plt.title("Mean forecast of BART model")
plt.show()
Your code seems to be syntactically correct, but I am not sure what the forecast variable represents. Essentially you are adding Gaussian noise to the mean function evaluated at the last observed value. Is that your intention?
Wow crazy 
You are really right. I do not know how it was possible that I did not see this! Thank you so much!
Will it be enough to put it like this and just specify mu with y?
# BART :-D
forecast_steps = 10
with pm.Model() as V_Research:
σ = pm.HalfNormal("σ", Y.std())
μ = pmb.BART("μ", X_lagged, Y, m=50)
y = pm.Normal("y", μ, σ, observed=Y)
forecast = pm.Normal("forecast", mu=y[-1], shape=(forecast_steps,))
idata_V_Research = pm.sample( return_inferencedata=True)
Now you are adding Gaussian noise to the mean predicted value evaluated at the last observed value. I don’t think that’s what you want. It seems you are looking for predictions (samples from the posterior predictive distribution). If that’s the case you have two options.
For in-sample predictions, you just need to call
pm.sample_posterior_predictive(idata_V_Research, V_Research, extend_inferencedata=True)
this will add a posterior_predictive group to idata_V_Research
For out-of-sample predictions, you instead need to define your model in a slightly different way.
with pm.Model() as V_Research:
Xs = pm.MutableData("Xs", X_lagged)
σ = pm.HalfNormal("σ", Y.std())
μ = pmb.BART("μ", Xs, Y, m=50)
y = pm.Normal("y", μ, σ, observed=Y, shape=μ.shape)
idata_V_Research = pm.sample()
and then
with V_Research:
Xs.set_value(X_test)
pm.sample_posterior_predictive(idata_V_Research, extend_inferencedata=True)
check Bayesian Additive Regression Trees: Introduction — PyMC-BART for a full example
You are right here as well. The mutable data is not the most easy thing for me.
Can I define X_test like below, and take the mean value of y for the direction of the forecast compared to the current value?
# Define X_test
X_test = np.zeros((10, X_lagged.shape[1]))
I am not sure I am following you. What do you mean by “for the direction of the forecast compared to the current value”?
Just to compare the Median of the forecast with the current value last value of the time series to get an idea of the future development.
I see, so you have something like X_train=[0, 1, 2, 3] and X_test=[4], and then you want to compare Y_train evaluated at X_train[3] and compare it with Y_test, right? If that’s the case notice that the value of Y estimated for X[3] will be the same for any value of X >=3.
The reason is that each leaf node returns a single value, then a sum of trees will be a piecewise constant function. Thus extrapolations (i.e. evaluation outside of the range of the training data) are going to be bad. This is true for BART and any other method that uses trees in a similar way.
pymc_bart.BART has a response argument that defaults to constant you can change it to linear, This will make each leaf node return a line. For some problems, this may help with extrapolations, for instance, if your unknown function is an increasing function, then the linear response should be able to capture that trend. If instead, your function is periodic, then the linear response will be of no help or it will work only for unobserved values close to the observed ones.
1 Like