def METHOD_NAME(self, n_samples, random_variables=None):
    """
    Execute the random sampling in the :class:`.KarhunenLoeveExpansion` class.
    The :meth:`run` method is the function that performs random sampling in the :class:`.KarhunenLoeveExpansion2D`
    class. If `n_samples` is provided when the :class:`.KarhunenLoeveExpansion2D` object is defined, the :meth:`run`
    method is automatically called. The user may also call the :meth:`run` method directly to generate samples.
    The :meth:`run`` method of the :class:`.KarhunenLoeveExpansion2D` class can be invoked many times and each time
    the generated samples are appended to the existing samples.
    :param n_samples: Number of samples of the stochastic process to be simulated.
     If the :meth:`run` method is invoked multiple times, the newly generated samples will be appended to the
     existing samples.
    The :meth:`run` method has no returns, although it creates and/or appends the :py:attr:`samples` attribute of
    the :class:`KarhunenLoeveExpansion2D` class.
    """
    samples = np.zeros((n_samples, self.correlation_function.shape[0], self.correlation_function.shape[2]))
    if random_variables is None:
        random_variables = np.random.normal(size=[self.thresholds[1], self.thresholds[0], n_samples])
    else:
        assert (random_variables.shape == (self.thresholds[1], self.thresholds[0], n_samples))
    for i in range(self.one_dimensional_correlation_function.shape[0]):
        if self.thresholds is not None:
            samples += np.einsum('x, xt, nx -> nxt', np.sqrt(self.w[:, i]), self.v[:, :, i],
                                 KarhunenLoeveExpansion(n_samples=n_samples,
                                                        correlation_function=
                                                        self.one_dimensional_correlation_function[i],
                                                        time_interval=self.time_intervals,
                                                        threshold=self.thresholds[0],
                                                        random_variables=random_variables[i]).samples[:, 0, :])
        else:
            samples += np.einsum('x, xt, nx -> nxt', np.sqrt(self.w[:, i]), self.v[:, :, i],
                                 KarhunenLoeveExpansion(n_samples=n_samples,
                                                        correlation_function=
                                                        self.one_dimensional_correlation_function[i],
                                                        time_interval=self.time_intervals,
                                                        random_variables=random_variables[i]).samples[:, 0, :])
    samples = np.reshape(samples, [samples.shape[0], 1, samples.shape[1], samples.shape[2]])
    if self.samples is None:
        self.samples = samples
        self.xi = random_variables
    else:
        self.samples = np.concatenate((self.samples, samples), axis=0)
        self.xi = np.concatenate((self.xi, random_variables), axis=2)