def sample(self, num_samples=1000, hmc_iters=20):
        """
        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
        """
        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