code: def predict_withGradients(self, X):
          
          if X.ndim==1: X = X[None,:]
          ps = self.model.param_array.copy()
          means = []
          stds = []
          dmdxs = []
          dsdxs = []
          for s in self.hmc_samples:
              if self.model._fixes_ is None:
                  self.model[:] = s
              else:
                  self.model[self.model._fixes_] = s
              self.model._trigger_params_changed()
              m, v = self.model.predict(X)
              std = np.sqrt(np.clip(v, 1e-10, np.inf))
              dmdx, dvdx = self.model.predictive_gradients(X)
              dmdx = dmdx[:,:,0]
              dsdx = dvdx / (2*std)
              means.append(m)
              stds.append(std)
              dmdxs.append(dmdx)
              dsdxs.append(dsdx)
          self.model.param_array[:] = ps
          self.model._trigger_params_changed()
          return means, stds, dmdxs, dsdxs
docs: Returns the mean, standard deviation, mean gradient and standard deviation gradient at X for all the MCMC samples.