code: def sample(self, num_samples=1000, hmc_iters=20):
          
          params = np.empty((num_samples,self.p.size))
          for i in range(num_samples):
              self.p[:] = np.random.multivariate_normal(np.zeros(self.p.size),self.M)
              H_old = self._computeH()
              theta_old = self.model.optimizer_array.copy()
              params[i] = self.model.unfixed_param_array
              #Matropolis
              self._update(hmc_iters)
              H_new = self._computeH()
  
              if H_old>H_new:
                  k = 1.
              else:
                  k = np.exp(H_old-H_new)
              if np.random.rand()<k:
                  params[i] = self.model.unfixed_param_array
              else:
                  self.model.optimizer_array = theta_old
          return params
docs: Sample the (unfixed) model parameters.
          
          :param num_samples: the number of samples to draw (1000 by default)
          :type num_samples: int
          :param hmc_iters: the number of leap-frog iterations (20 by default)
          :type hmc_iters: int
          :return: the list of parameters samples with the size N x P (N - the number of samples, P - the number of parameters to sample) 
          :rtype: numpy.ndarray