It seems I am not to far from a result:
# Create the model
with pm.Model() as model_disasters:
μ_ = pmb.BART("μ_", X=x_data, Y=y_data, m=20)
μ = pm.Deterministic("μ", pm.math.abs(μ_))
y_pred = pm.NegativeBinomial("y_pred", mu=μ, alpha=1, observed=y_data)
idata_disasters = pm.sample(random_seed=RANDOM_SEED)
# Forecast 10 steps into the future
future_years = 10
x_future = np.arange(x_data.min(), x_data.max() + future_years + 1)[:, None]
y_future = np.zeros((x_future.shape[0],))
with model_disasters:
# Generate samples from the posterior predictive distribution
post_pred = pm.sample_posterior_predictive(idata_disasters, var_names=["y_pred"], random_seed=RANDOM_SEED)
# Take the mean of the posterior predictive distribution as the forecast
y_future = post_pred["y_pred"].mean(axis=0)
# Plot the forecast
_, ax = plt.subplots(figsize=(10, 6))
ax.plot(x_centers, disaster_rate_mean, "w", lw=3)
az.plot_hdi(x_centers, disaster_rates, smooth=False)
az.plot_hdi(x_centers, disaster_rates, hdi_prob=0.5, smooth=False, plot_kwargs={"alpha": 0})
ax.plot(disasters, np.zeros_like(disasters) - 0.5, "k|")
ax.plot(x_future, y_future, 'g--', label='forecast')
ax.set_xlabel("years")
ax.set_ylabel("disaster_rate")
ax.legend()
plt.show()
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Input In [16], in <cell line: 13>()
11 y_future = np.zeros((x_future.shape[0],))
13 with model_disasters:
14 # Generate samples from the posterior predictive distribution
---> 15 post_pred = pm.sample_posterior_predictive(idata_disasters, var_names=["y_pred"], random_seed=RANDOM_SEED)
16 # Take the mean of the posterior predictive distribution as the forecast
17 y_future = post_pred["y_pred"].mean(axis=0)
File ~\anaconda3\lib\site-packages\pymc\sampling\forward.py:673, in sample_posterior_predictive(trace, model, var_names, sample_dims, random_seed, progressbar, return_inferencedata, extend_inferencedata, predictions, idata_kwargs, compile_kwargs)
671 ikwargs.setdefault("inplace", True)
672 return pm.predictions_to_inference_data(ppc_trace, **ikwargs)
--> 673 idata_pp = pm.to_inference_data(posterior_predictive=ppc_trace, **ikwargs)
675 if extend_inferencedata and idata is not None:
676 idata.extend(idata_pp)
File ~\anaconda3\lib\site-packages\pymc\backends\arviz.py:485, in to_inference_data(trace, prior, posterior_predictive, log_likelihood, coords, dims, sample_dims, model, save_warmup, include_transformed)
482 if isinstance(trace, InferenceData):
483 return trace
--> 485 return InferenceDataConverter(
486 trace=trace,
487 prior=prior,
488 posterior_predictive=posterior_predictive,
489 log_likelihood=log_likelihood,
490 coords=coords,
491 dims=dims,
492 sample_dims=sample_dims,
493 model=model,
494 save_warmup=save_warmup,
495 include_transformed=include_transformed,
496 ).to_inference_data()
File ~\anaconda3\lib\site-packages\pymc\backends\arviz.py:406, in InferenceDataConverter.to_inference_data(self)
396 def to_inference_data(self):
397 """Convert all available data to an InferenceData object.
398
399 Note that if groups can not be created (e.g., there is no `trace`, so
400 the `posterior` and `sample_stats` can not be extracted), then the InferenceData
401 will not have those groups.
402 """
403 id_dict = {
404 "posterior": self.posterior_to_xarray(),
405 "sample_stats": self.sample_stats_to_xarray(),
--> 406 "posterior_predictive": self.posterior_predictive_to_xarray(),
407 "predictions": self.predictions_to_xarray(),
408 **self.priors_to_xarray(),
409 "observed_data": self.observed_data_to_xarray(),
410 }
411 if self.predictions:
412 id_dict["predictions_constant_data"] = self.constant_data_to_xarray()
File ~\anaconda3\lib\site-packages\arviz\data\base.py:65, in requires.__call__.<locals>.wrapped(cls)
63 if all((getattr(cls, prop_i) is None for prop_i in prop)):
64 return None
---> 65 return func(cls)
File ~\anaconda3\lib\site-packages\pymc\backends\arviz.py:327, in InferenceDataConverter.posterior_predictive_to_xarray(self)
325 data = self.posterior_predictive
326 dims = {var_name: self.sample_dims + self.dims.get(var_name, []) for var_name in data}
--> 327 return dict_to_dataset(
328 data, library=pymc, coords=self.coords, dims=dims, default_dims=self.sample_dims
329 )
File ~\anaconda3\lib\site-packages\arviz\data\base.py:307, in dict_to_dataset(data, attrs, library, coords, dims, default_dims, index_origin, skip_event_dims)
305 data_vars = {}
306 for key, values in data.items():
--> 307 data_vars[key] = numpy_to_data_array(
308 values,
309 var_name=key,
310 coords=coords,
311 dims=dims.get(key),
312 default_dims=default_dims,
313 index_origin=index_origin,
314 skip_event_dims=skip_event_dims,
315 )
316 return xr.Dataset(data_vars=data_vars, attrs=make_attrs(attrs=attrs, library=library))
File ~\anaconda3\lib\site-packages\arviz\data\base.py:254, in numpy_to_data_array(ary, var_name, coords, dims, default_dims, index_origin, skip_event_dims)
252 # filter coords based on the dims
253 coords = {key: xr.IndexVariable((key,), data=np.asarray(coords[key])) for key in dims}
--> 254 return xr.DataArray(ary, coords=coords, dims=dims)
File ~\anaconda3\lib\site-packages\xarray\core\dataarray.py:428, in DataArray.__init__(self, data, coords, dims, name, attrs, indexes, fastpath)
426 data = _check_data_shape(data, coords, dims)
427 data = as_compatible_data(data)
--> 428 coords, dims = _infer_coords_and_dims(data.shape, coords, dims)
429 variable = Variable(dims, data, attrs, fastpath=True)
430 indexes, coords = _create_indexes_from_coords(coords)
File ~\anaconda3\lib\site-packages\xarray\core\dataarray.py:142, in _infer_coords_and_dims(shape, coords, dims)
140 dims = tuple(dims)
141 elif len(dims) != len(shape):
--> 142 raise ValueError(
143 "different number of dimensions on data "
144 f"and dims: {len(shape)} vs {len(dims)}"
145 )
146 else:
147 for d in dims:
ValueError: different number of dimensions on data and dims: 3 vs 2
If I find out how to put the mutable data the model should be ready and running without errors. I am using Chat GPT as an assistent, but it is only trained up to data from 2021. So not on the new version of BART v. 0.4.0 which still makes it tricky for me 
but for today I give up …
You got it right that I try to reproduce the disscusion and the advise from Osvaldo, but in his solution there is nothing about obs…
