def draw_sample(self, Xstar, n=0, num_samp=1, rand_vars=None,
                    rand_type='standard normal', diag_factor=1e3,
                    method='cholesky', num_eig=None, mean=None, cov=None,
                    modify_sign=None, **kwargs):
        """Draw a sample evaluated at the given points `Xstar`.
        
        Note that this function draws samples from the GP given the current
        values for the hyperparameters (which may be in a nonsense state if you
        just created the instance or called a method that performs MCMC sampling).
        If you want to draw random samples from MCMC output, use the
        `return_samples` and `full_output` keywords to :py:meth:`predict`.
        
        Parameters
        ----------
        Xstar : array, (`M`, `D`)
            `M` test input values of dimension `D`.
        n : array, (`M`, `D`) or scalar, non-negative int, optional
            Derivative order to evaluate at. Default is 0 (evaluate value).
        noise : bool, optional
            Whether or not to include the noise components of the kernel in the
            sample. Default is False (no noise in samples).
        num_samp : Positive int, optional
            Number of samples to draw. Default is 1. Cannot be used in
            conjunction with `rand_vars`: If you pass both `num_samp` and
            `rand_vars`, `num_samp` will be silently ignored.
        rand_vars : array, (`M`, `P`), optional
            Vector of random variables :math:`u` to use in constructing the
            sample :math:`y_* = f_* + Lu`, where :math:`K=LL^T`. If None,
            values will be produced using
            :py:func:`numpy.random.multivariate_normal`. This allows you to use
            pseudo/quasi random numbers generated by an external routine. Note
            that, when `method` is 'eig', the eigenvalues are in *ascending*
            order.
            Default is None (use :py:func:`multivariate_normal` directly).
        rand_type : {'standard normal', 'uniform'}, optional
            Type of distribution the inputs are given with.
            
                * 'standard normal': Standard (`mu` = 0, `sigma` = 1) normal
                  distribution (this is the default)
                * 'uniform': Uniform distribution on [0, 1). In this case
                  the required Gaussian variables are produced with inversion.
                  
        diag_factor : float, optional
            Number (times machine epsilon) added to the diagonal of the
            covariance matrix prior to computing its Cholesky decomposition.
            This is necessary as sometimes the decomposition will fail because,
            to machine precision, the matrix appears to not be positive definite.
            If you are getting errors from :py:func:`scipy.linalg.cholesky`, try
            increasing this an order of magnitude at a time. This parameter only
            has an effect when using rand_vars. Default value is 1e3. 
        method : {'cholesky', 'eig'}, optional
            Method to use for constructing the matrix square root. Default is
            'cholesky' (use lower-triangular Cholesky decomposition).
            
                * 'cholesky': Perform Cholesky decomposition on the covariance
                  matrix: :math:`K=LL^T`, use :math:`L` as the matrix square
                  root.
                * 'eig': Perform an eigenvalue decomposition on the covariance
                  matrix: :math:`K=Q \\Lambda Q^{-1}`, use :math:`Q\\Lambda^{1/2}`
                  as the matrix square root.
        
        num_eig : int or None, optional
            Number of eigenvalues to compute. Can range from 1 to `M` (the
            number of test points). If it is None, then all eigenvalues are
            computed. Default is None (compute all eigenvalues). This keyword
            only has an effect if `method` is 'eig'.
        mean : array, (`M`,), optional
            If you have pre-computed the mean and covariance matrix, then you
            can simply pass them in with the `mean` and `cov` keywords to save
            on having to call :py:meth:`predict`.
        cov : array, (`M`, `M`), optional
            If you have pre-computed the mean and covariance matrix, then you
            can simply pass them in with the `mean` and `cov` keywords to save
            on having to call :py:meth:`predict`.
        modify_sign : {None, 'left value', 'right value', 'left slope', 'right slope', 'left concavity', 'right concavity'}, optional
            If None (the default), the eigenvectors as returned by
            :py:func:`scipy.linalg.eigh` are used without modification. To
            modify the sign of the eigenvectors (necessary for some advanced use
            cases), set this kwarg to one of the following:
            
                * 'left value': forces the first value of each eigenvector to be
                  positive.
                * 'right value': forces the last value of each eigenvector to be
                  positive.
                * 'left slope': forces the slope to be positive at the start of
                  each eigenvector.
                * 'right slope': forces the slope to be positive at the end of
                  each eigenvector.
                * 'left concavity': forces the second derivative to be positive
                  at the start of each eigenvector.
                * 'right concavity': forces the second derivative to be positive
                  at the end of each eigenvector.
        
        **kwargs : optional kwargs
            All extra keyword arguments are passed to :py:meth:`predict` when
            evaluating the mean and covariance matrix of the GP.
        
        Returns
        -------
        samples : :py:class:`Array` (`M`, `P`) or (`M`, `num_samp`)
            Samples evaluated at the `M` points.
        
        Raises
        ------
        ValueError
            If rand_type or method is invalid.
        """
        # All of the input processing for Xstar and n will be done in here:
        if mean is None or cov is None:
            out = self.predict(Xstar, n=n, full_output=True, **kwargs)
            mean = out['mean']
            cov = out['cov']
        if rand_vars is None and method != 'eig':
            try:
                return numpy.random.multivariate_normal(mean, cov, num_samp).T
            except numpy.linalg.LinAlgError as e:
                if self.verbose:
                    warnings.warn(
                        "Failure when drawing from MVN! Falling back on eig. "
                        "Exception was:\n%s"
                        % (e,),
                        RuntimeWarning
                    )
                method = 'eig'
        
        if num_eig is None or num_eig > len(mean):
            num_eig = len(mean)
        elif num_eig < 1:
            num_eig = 1
        if rand_vars is None:
            rand_vars = numpy.random.standard_normal((num_eig, num_samp))
        valid_types = ('standard normal', 'uniform')
        if rand_type not in valid_types:
            raise ValueError(
                "rand_type %s not recognized! Valid options are: %s."
                % (rand_type, valid_types,)
            )
        if rand_type == 'uniform':
            rand_vars = scipy.stats.norm.ppf(rand_vars)
        
        if method == 'cholesky':
            L = scipy.linalg.cholesky(
                cov + diag_factor * sys.float_info.epsilon * scipy.eye(cov.shape[0]),
                lower=True,
                check_finite=False
            )
        elif method == 'eig':
            # TODO: Add support for specifying cutoff eigenvalue!
            # Not technically lower triangular, but we'll keep the name L:
            eig, Q = scipy.linalg.eigh(
                cov + diag_factor * sys.float_info.epsilon * scipy.eye(cov.shape[0]),
                eigvals=(len(mean) - 1 - (num_eig - 1), len(mean) - 1)
            )
            if modify_sign is not None:
                if modify_sign == 'left value':
                    modify_mask = (Q[0, :] < 0.0)
                elif modify_sign == 'right value':
                    modify_mask = (Q[-1, :] < 0.0)
                elif modify_sign == 'left slope':
                    modify_mask = ((Q[1, :] - Q[0, :]) < 0.0)
                elif modify_sign == 'right slope':
                    modify_mask = ((Q[-1, :] - Q[-2, :]) < 0.0)
                elif modify_sign == 'left concavity':
                    modify_mask = ((Q[2, :] - 2 * Q[1, :] + Q[0, :]) < 0.0)
                elif modify_sign == 'right concavity':
                    modify_mask = ((Q[-1, :] - 2 * Q[-2, :] + Q[-3, :]) < 0.0)
                else:
                    raise ValueError(
                        "modify_sign %s not recognized!" % (modify_sign,)
                    )
                Q[:, modify_mask] *= -1.0
            Lam_1_2 = scipy.diag(scipy.sqrt(eig))
            L = Q.dot(Lam_1_2)
        else:
            raise ValueError("method %s not recognized!" % (method,))
        return scipy.atleast_2d(mean).T + L.dot(rand_vars[:num_eig, :])