I received an error when running the following code.
the error is “input operand has more dimensions than allowed by the axis remapping”

# set the seed
np.random.seed(1)
n = 100 # The number of data points
X = np.linspace(0, 10, n)[:, None] # The inputs to the GP, they must be arranged as a column vector
# Define the true covariance function and its parameters
ℓ_true = 1.0
η_true = 3.0
cov_func = η_true**2 * pm.gp.cov.Matern52(1, ℓ_true)
# A mean function that is zero everywhere
mean_func = pm.gp.mean.Zero()
# The latent function values are one sample from a multivariate normal
# Note that we have to call `eval()` because PyMC3 built on top of Theano
tp_samples = pm.MvStudentT.dist(mu=mean_func(X), cov=cov_func(X), nu=3).random(size=8)
## Plot samples from TP prior
fig = plt.figure(figsize=(12, 5))
ax = fig.gca()
ax.plot(X.flatten(), tp_samples.T, lw=3, alpha=0.6)
ax.set_xlabel("X")
ax.set_ylabel("y")
ax.set_title("Samples from TP with DoF=3")
gp_samples = pm.MvNormal.dist(mu=mean_func(X).eval(), cov=cov_func(X).eval()).random(size=8)
fig = plt.figure(figsize=(12, 5))
ax = fig.gca()
ax.plot(X.flatten(), gp_samples.T, lw=3, alpha=0.6)
ax.set_xlabel("X")
ax.set_ylabel("y")
ax.set_title("Samples from GP");

Either your random size is invalid, or there is a bug. This would very likely work out of the box in pymc>4.0, so I keep my suggestion of trying to update it.