GP-MeansAndCovs.ipynb throws LinAlgError

Cell 6 contains:

lengthscale = 0.2
eta = 2.0
cov = eta ** 2 * pm.gp.cov.ExpQuad(1, lengthscale)

X = np.linspace(0, 2, 200)[:, None]
K = cov(X).eval()

plt.figure(figsize=(14, 4))

plt.plot(X, pm.MvNormal.dist(mu=np.zeros(K.shape[0]), cov=K).random(size=3).T)
plt.title("Samples from the GP prior")
plt.ylabel("y")
plt.xlabel("X");

which throws

ValueError: input operand has more dimensions than allowed by the axis remapping

As per #11 , I can fix this by specifying the shape of the MvNormal:

lengthscale = 0.2
eta = 2.0
cov = eta ** 2 * pm.gp.cov.ExpQuad(1, lengthscale)

X = np.linspace(0, 2, 200)[:, None]
K = cov(X).eval()

plt.figure(figsize=(14, 4))

plt.plot(X, pm.MvNormal.dist(mu=np.zeros(K.shape[0]), cov=K, shape=K.shape[0]).random(size=3).T)
plt.title("Samples from the GP prior")
plt.ylabel("y")
plt.xlabel("X");

but then I get

~/pymc3-dev/pymc3/distributions/multivariate.py in random(self, point, size)
    277 
    278         if self._cov_type == "cov":
--> 279             chol = np.linalg.cholesky(param)
    280         elif self._cov_type == "chol":
    281             chol = param

<__array_function__ internals> in cholesky(*args, **kwargs)

~/miniconda3/envs/pymc3-dev-py38/lib/python3.8/site-packages/numpy/linalg/linalg.py in cholesky(a)
    762     t, result_t = _commonType(a)
    763     signature = 'D->D' if isComplexType(t) else 'd->d'
--> 764     r = gufunc(a, signature=signature, extobj=extobj)
    765     return wrap(r.astype(result_t, copy=False))
    766 

~/miniconda3/envs/pymc3-dev-py38/lib/python3.8/site-packages/numpy/linalg/linalg.py in _raise_linalgerror_nonposdef(err, flag)
     89 
     90 def _raise_linalgerror_nonposdef(err, flag):
---> 91     raise LinAlgError("Matrix is not positive definite")
     92 
     93 def _raise_linalgerror_eigenvalues_nonconvergence(err, flag):

LinAlgError: Matrix is not positive definite

cc @Sayam753 any suggestions?

Open GP-MeansAndCovs.ipynb and reproduce Cell 6 in the reported environment, first without and then with the MvNormal shape argument. Trace the covariance construction and sampling path around the reported errors; done means the cell runs successfully and produces the GP prior sample plot without either exception.