Hi, I’m trying to calibrate a PyMC Bayesian model with new data. Basically I would like to update the prior with the posterior distributions of the pre-trained model. I followed a similar workflow of this work: Updating priors — PyMC3 3.11.5 documentation. I understand general idea is to find the gaussian kernel density estimation of the posterior samples and rebuild corresponding priors with `pymc.Interpolated`

distribution.

The only difference between my work and the demo notebook is that some of the coefficients in my model are multi-dimensional. I parameterized the argument `shape`

and pass an array to distribution parameters, such as `mu`

and `sigma`

, to build a multi-dimensional `Normal`

distribution object. See example codes below:

```
import pymc3 as pm
import numpy as np
import theano.tensor as tt
model = pm.Model()
with model:
...
beta = pm.Normal(
name="beta",
mu=tt._shared(np.array([1, 2, 3])),
sigma=tt._shared(np.array([1, 2, 3])),
shape=3
)
...
trace = pm.sample(draws=2000, tune=2000, ...)
```

Then I modified the `from_posterior`

function with `np.apply_along_axis`

to generate a batch of `x_points`

and `pdf_points`

.So that I can similarly pass arrays of shape `(3, 1000)`

to argument `x_points`

and `pdf_points`

, where the slice along axis 1 parameterize each dimension of the distribution. See code examples below:

```
import pymc3 as pm
from scipy import stats
import numpy as np
# I modified the function a little bit because my sampled posteriors are multi-dimensional
def from_posterior(tag, samples, bins=1000):
if len(samples.shape) == 1:
samples = np.expand_dims(samples, axis=1)
smin, smax = np.min(samples, axis=0), np.max(samples, axis=0)
width = smax - smin
x = np.linspace(smin, smax, bins)
y= np.apply_along_axis(func1d=lambda a: stats.gaussian_kde(a[:-bins])(a[-bins:]),
axis=0,
arr=np.concatenate((samples, x), axis=0))
x_points = np.concatenate((x[[0], :] - 3 * width, x, x[[-1], :] + 3 * width), axis=0).T
# x_points.shape = (samples.shape[1], bins+2)
pdf_points = np.concatenate([np.zeros((1, y.shape[1])), y, np.zeros((1, y.shape[1]))], axis=0).T
# pdf_points.shape = (samples.shape[1], bins+2)
return pm.Interpolated(tag, x_points=x_points, pdf_points=pdf_points, shape=samples.shape[1])
with model:
...
beta = from_posterior(
tag="beta",
samples=beta_posteriors, # beta_posteriors.shape = (8000, 3)
bins=1000
)
...
```

However, the codes above returned this error:

```
ValueError: too many axes: 2 (effrank=2), expected rank=1
```

It seems the class `pymc.Interpolated`

doesn’t function similarly as `pm.Normal`

. The source codes called `scipy.interpolate`

which doesn’t support the syntax above. So I just wonder if anyone encountered a similar problem and knows a way to code the dimension size of `pm.Interpolated`

dynamically? Thank you very much!