diff --git a/docs/code/surrogates/gpr/plot_gpr_custom2D.py b/docs/code/surrogates/gpr/plot_gpr_custom2D.py
index 97820ea8c..764130894 100644
--- a/docs/code/surrogates/gpr/plot_gpr_custom2D.py
+++ b/docs/code/surrogates/gpr/plot_gpr_custom2D.py
@@ -112,7 +112,7 @@
 y_act = np.array(r2model.qoi_list).reshape(x1g.shape[0], x1g.shape[1])
 
 fig1 = plt.figure()
-ax = fig1.gca(projection='3d')
+ax = fig1.add_subplot(projection='3d')
 surf = ax.plot_surface(x1g, x2g, y_act, cmap=cm.coolwarm, linewidth=0, antialiased=False)
 ax.set_zlim(-1, 15)
 ax.zaxis.set_major_locator(LinearLocator(10))
diff --git a/docs/code/surrogates/pce/plot_pce_camel.py b/docs/code/surrogates/pce/plot_pce_camel.py
index 27c3e7ed6..7e0a0cafb 100644
--- a/docs/code/surrogates/pce/plot_pce_camel.py
+++ b/docs/code/surrogates/pce/plot_pce_camel.py
@@ -76,7 +76,7 @@ def function(x, y):
 f = function(X1_, X2_)
 
 fig = plt.figure(figsize=(10, 6))
-ax = fig.gca(projection='3d')
+ax = fig.add_subplot(projection='3d')
 surf = ax.plot_surface(X1_, X2_, f, rstride=1, cstride=1, cmap='gnuplot2', linewidth=0, antialiased=False)
 ax.set_title('True function')
 ax.set_xlabel('$x_1$', fontsize=15)
@@ -95,7 +95,7 @@ def function(x, y):
 # %%
 
 fig = plt.figure(figsize=(10, 6))
-ax = fig.gca(projection='3d')
+ax = fig.add_subplot(projection='3d')
 ax.scatter(x[:, 0], x[:, 1], y, s=20, c='r')
 
 ax.set_title('Training data')
@@ -168,7 +168,7 @@ def function(x, y):
 # %%
 
 fig = plt.figure(figsize=(10,6))
-ax = fig.gca(projection='3d')
+ax = fig.add_subplot(projection='3d')
 ax.scatter(x_test[:,0], x_test[:,1], y_test, s=1)
 
 ax.set_title('PCE predictor')
diff --git a/docs/code/surrogates/pce/plot_pce_sphere.py b/docs/code/surrogates/pce/plot_pce_sphere.py
index 57382c7e9..1c6e42dde 100644
--- a/docs/code/surrogates/pce/plot_pce_sphere.py
+++ b/docs/code/surrogates/pce/plot_pce_sphere.py
@@ -71,7 +71,7 @@ def function(x,y):
 f = function(X1_, X2_)
 
 fig = plt.figure(figsize=(10,6))
-ax = fig.gca(projection='3d')
+ax = fig.add_subplot(projection='3d')
 surf = ax.plot_surface(X1_, X2_, f, rstride=1, cstride=1, cmap='gnuplot2', linewidth=0, antialiased=False)
 ax.set_title('True function')
 ax.set_xlabel('$x_1$', fontsize=15)
@@ -90,7 +90,7 @@ def function(x,y):
 # %%
 
 fig = plt.figure(figsize=(10,6))
-ax = fig.gca(projection='3d')
+ax = fig.add_subplot(projection='3d')
 ax.scatter(x[:,0], x[:,1], y, s=20, c='r')
 
 ax.set_title('Training data')
@@ -156,7 +156,7 @@ def function(x,y):
 # %%
 
 fig = plt.figure(figsize=(10,6))
-ax = fig.gca(projection='3d')
+ax = fig.add_subplot(projection='3d')
 ax.scatter(x_test[:,0], x_test[:,1], y_test, s=1)
 
 ax.set_title('PCE predictor')
diff --git a/src/UQpy/sampling/adaptive_kriging_functions/baseclass/LearningFunction.py b/src/UQpy/sampling/adaptive_kriging_functions/baseclass/LearningFunction.py
index 1da00c410..e754dfb01 100644
--- a/src/UQpy/sampling/adaptive_kriging_functions/baseclass/LearningFunction.py
+++ b/src/UQpy/sampling/adaptive_kriging_functions/baseclass/LearningFunction.py
@@ -4,7 +4,7 @@
 class LearningFunction(ABC):
     def __init(self, ordered_parameters=None, **kwargs):
         self.parameters = kwargs
-        self.ordered_parameters = (ordered_parameters if not None else tuple(kwargs.keys()))
+        self.ordered_parameters = (ordered_parameters if ordered_parameters is not None else tuple(kwargs.keys()))
         if len(self.ordered_parameters) != len(self.parameters):
             raise ValueError("Inconsistent dimensions between order_params tuple and params dictionary.")
 
diff --git a/src/UQpy/transformations/Nataf.py b/src/UQpy/transformations/Nataf.py
index 54c4075ca..90c3adf2c 100644
--- a/src/UQpy/transformations/Nataf.py
+++ b/src/UQpy/transformations/Nataf.py
@@ -21,17 +21,17 @@ class Nataf:
 
     @beartype
     def __init__(
-        self,
-        distributions: Union[Distribution, DistributionList],
-        samples_x: Union[None, np.ndarray] = None,
-        samples_z: Union[None, np.ndarray] = None,
-        jacobian: bool = False,
-        corr_z: Union[None, np.ndarray] = None,
-        corr_x: Union[None, np.ndarray] = None,
-        itam_beta: Union[float, int] = 1.0,
-        itam_threshold1: Union[float, int] = 0.001,
-        itam_threshold2: Union[float, int] = 0.1,
-        itam_max_iter: int = 100,
+            self,
+            distributions: Union[Distribution, DistributionList],
+            samples_x: Union[None, np.ndarray] = None,
+            samples_z: Union[None, np.ndarray] = None,
+            jacobian: bool = False,
+            corr_z: Union[None, np.ndarray] = None,
+            corr_x: Union[None, np.ndarray] = None,
+            itam_beta: Union[float, int] = 1.0,
+            itam_threshold1: Union[float, int] = 0.001,
+            itam_threshold2: Union[float, int] = 0.1,
+            itam_max_iter: int = 100,
     ):
         """
         Transform random variables using the Nataf or Inverse Nataf transformation
@@ -73,7 +73,7 @@ def __init__(
         self.dist_object = distributions
         self.samples_x: NumpyFloatArray = samples_x
         """Random vector of shape ``(nsamples, n_dimensions)`` with prescribed probability distributions."""
-        self.samples_z:NumpyFloatArray = samples_z
+        self.samples_z: NumpyFloatArray = samples_z
         """Standard normal random vector of shape ``(nsamples, n_dimensions)``"""
         self.jacobian = jacobian
         self.jzx: NumpyFloatArray = None
@@ -98,9 +98,9 @@ def __init__(
             elif all(isinstance(x, Normal) for x in distributions):
                 self.corr_z = self.corr_x
             else:
-                self.corr_z, self.itam_error1, self.itam_error2 =\
+                self.corr_z, self.itam_error1, self.itam_error2 = \
                     self.itam(self.dist_object, self.corr_x, self.itam_max_iter, self.itam_beta,
-                              self.itam_threshold1, self.itam_threshold2,)
+                              self.itam_threshold1, self.itam_threshold2, )
         elif corr_z is not None:
             self.corr_z = corr_z
             if np.all(np.equal(self.corr_z, np.eye(self.n_dimensions))):
@@ -108,7 +108,7 @@ def __init__(
             elif all(isinstance(x, Normal) for x in distributions):
                 self.corr_x = self.corr_z
             else:
-                self.corr_x = self.distortion_z2x(self.dist_object, self.corr_z)
+                self.corr_x = self.distortion_z2x(self.dist_object, self.corr_z, n_gauss_points=128)
 
         self.H: NumpyFloatArray = cholesky(self.corr_z, lower=True)
         """The lower triangular matrix resulting from the Cholesky decomposition of the correlation matrix
@@ -119,10 +119,10 @@ def __init__(
 
     @beartype
     def run(
-        self,
-        samples_x: Union[None, np.ndarray] = None,
-        samples_z: Union[None, np.ndarray] = None,
-        jacobian: bool = False,
+            self,
+            samples_x: Union[None, np.ndarray] = None,
+            samples_z: Union[None, np.ndarray] = None,
+            jacobian: bool = False,
     ):
         """
         Execute the Nataf transformation or its inverse.
@@ -160,15 +160,15 @@ def run(
 
     @staticmethod
     def itam(
-        distributions: Union[
-            DistributionContinuous1D,
-            JointIndependent,
-            list[Union[DistributionContinuous1D, JointIndependent]]],
-        corr_x,
-        itam_max_iter: int = 100,
-        itam_beta: Union[float, int] = 1.0,
-        itam_threshold1: Union[float, int] = 0.001,
-        itam_threshold2: Union[float, int] = 0.01,
+            distributions: Union[
+                DistributionContinuous1D,
+                JointIndependent,
+                list[Union[DistributionContinuous1D, JointIndependent]]],
+            corr_x,
+            itam_max_iter: int = 100,
+            itam_beta: Union[float, int] = 1.0,
+            itam_threshold1: Union[float, int] = 0.001,
+            itam_threshold2: Union[float, int] = 0.01,
     ):
         """
         Calculate the correlation matrix :math:`\mathbf{C_Z}` of the standard normal random vector
@@ -236,7 +236,8 @@ def itam(
         return corr_z, itam_error1, itam_error2
 
     @staticmethod
-    def distortion_z2x(distributions:  Union[Distribution, list[Distribution]], corr_z: np.ndarray):
+    def distortion_z2x(distributions: Union[Distribution, list[Distribution]], corr_z: np.ndarray,
+                       n_gauss_points: int = 1024):
         """
         This is a method to calculate the correlation matrix :math:`\mathbf{C_x}` of the random vector
         :math:`\mathbf{x}`  given the correlation matrix :math:`\mathbf{C_z}` of the standard normal random vector
@@ -248,12 +249,14 @@ def distortion_z2x(distributions:  Union[Distribution, list[Distribution]], corr
          This method is part of the :class:`.Nataf` class.
         :param corr_z: The correlation  matrix (:math:`\mathbf{C_z}`) of the standard normal vector **Z** .
          Default: The ``identity`` matrix.
+        :param n_gauss_points: The number of integration points used for the numerical integration of the
+         correlation  matrix (:math:`\mathbf{C_Z}`) of the standard normal random vector **Z**
         :return: Distorted correlation matrix (:math:`\mathbf{C_X}`) of the random vector **X**.
         """
         logger = logging.getLogger(__name__)
         z_max = 8
         z_min = -z_max
-        ng = 128
+        ng = n_gauss_points
 
         eta, w2d, xi = calculate_gauss_quadrature_2d(ng, z_max, z_min)
 
@@ -268,7 +271,7 @@ def distortion_z2x(distributions:  Union[Distribution, list[Distribution]], corr
     @staticmethod
     def calculate_corr_x(corr_x, corr_z, marginals, eta, w2d, xi, is_joint):
         if all(hasattr(m, "moments") for m in marginals) and all(
-            hasattr(m, "icdf") for m in marginals
+                hasattr(m, "icdf") for m in marginals
         ):
             for i in range(len(marginals)):
                 i_cdf_i = marginals[i].icdf