diff --git a/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu_rqmc_based.py b/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu_rqmc_based.py
index 2f3da81..fc089cf 100644
--- a/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu_rqmc_based.py
+++ b/src/algorithms/semiparam_algorithms/nvm_semi_param_algorithms/g_estimation_given_mu_rqmc_based.py
@@ -7,7 +7,8 @@
 from scipy.integrate import quad_vec
 from scipy.special import gamma
 
-from src.algorithms.support_algorithms.rqmc import RQMC
+from src.algorithms.support_algorithms.integrator import Integrator
+from src.algorithms.support_algorithms.rqmc import RQMCIntegrator
 from src.estimators.estimate_result import EstimateResult
 
 MU_DEFAULT_VALUE = 1.0
@@ -68,7 +69,7 @@ def __init__(self, sample: Optional[_typing.NDArray[np.float64]] = None, **kwarg
             self.x_data,
             self.grid_size,
             self.integration_tolerance,
-            self.integration_limit,
+            self.integration_limit
         ) = self._validate_kwargs(self.n, **kwargs)
         self.denominator: float = 2 * math.pi * self.n
         self.precompute_gamma_grid()
@@ -162,15 +163,15 @@ def second_v_integrand(self, v: float, x: float) -> np.ndarray:
         x_power = self.x_powers[x][idx]
         return (self.second_u_integrals[idx] * x_power) / gamma_val
 
-    def compute_integrals_for_x(self, x: float) -> float:
+    def compute_integrals_for_x(self, x: float, integrator: Integrator = RQMCIntegrator()) -> float:
         """Compute integrals using RQMC for v-integration."""
-        first_integral = RQMC(lambda t: np.sum(self.first_v_integrand(t * self.v_value, x)) * self.v_value).rqmc()[0]
+        first_integral = integrator.compute(func=lambda t: np.sum(self.first_v_integrand(t * self.v_value, x)) * self.v_value).value
 
-        second_integral = RQMC(lambda t: np.sum(self.second_v_integrand(-t * self.v_value, x)) * self.v_value).rqmc()[0]
+        second_integral = integrator.compute(func=lambda t: np.sum(self.second_v_integrand(-t * self.v_value, x)) * self.v_value).value
 
         total = (first_integral + second_integral) / self.denominator
         return max(0.0, total.real)
 
     def algorithm(self, sample: np._typing.NDArray) -> EstimateResult:
-        y_data = [self.compute_integrals_for_x(x) for x in self.x_data]
+        y_data = [self.compute_integrals_for_x(x, RQMCIntegrator()) for x in self.x_data]
         return EstimateResult(list_value=y_data, success=True)
diff --git a/src/algorithms/support_algorithms/integrator.py b/src/algorithms/support_algorithms/integrator.py
new file mode 100644
index 0000000..8ee0fb7
--- /dev/null
+++ b/src/algorithms/support_algorithms/integrator.py
@@ -0,0 +1,20 @@
+from dataclasses import dataclass
+from typing import Any, Protocol, Callable, Optional
+
+
+@dataclass
+class IntegrationResult:
+    value: float
+    error: float
+    message: Optional[dict[str, Any]] | None = None
+
+
+class Integrator(Protocol):
+
+    """Base class for integral calculation"""
+
+    def __init__(self) -> None:
+        ...
+
+    def compute(self, func: Callable) -> IntegrationResult:
+        ...
diff --git a/src/algorithms/support_algorithms/quad_integrator.py b/src/algorithms/support_algorithms/quad_integrator.py
new file mode 100644
index 0000000..616eca2
--- /dev/null
+++ b/src/algorithms/support_algorithms/quad_integrator.py
@@ -0,0 +1,60 @@
+from typing import Callable, Any, Sequence
+
+from scipy.integrate import quad
+from src.algorithms.support_algorithms.integrator import IntegrationResult
+
+class QuadIntegrator:
+
+    def __init__(
+            self,
+            a: float = 0,
+            b: float = 1,
+            args: tuple[Any, ...] = (),
+            full_output: int = 0,
+            epsabs: float | int = 1.49e-08,
+            epsrel: float | int = 1.49e-08,
+            limit: float | int = 50,
+            points: Sequence[float | int] | None = None,
+            weight: float | int | None = None,
+            wvar: Any = None,
+            wopts: Any = None,
+            maxp1: float | int = 50,
+            limlst: int = 50,
+            complex_func: bool = False,
+    ):
+        self.params = {
+            'a': a,
+            'b': b,
+            'args': args,
+            'full_output': full_output,
+            'epsabs': epsabs,
+            'epsrel': epsrel,
+            'limit': limit,
+            'points': points,
+            'weight': weight,
+            'wvar': wvar,
+            'wopts': wopts,
+            'maxp1': maxp1,
+            'limlst': limlst,
+            'complex_func': complex_func
+        }
+
+    def compute(self, func: Callable) -> IntegrationResult:
+
+        """
+        Compute integral via quad integrator
+
+        Args:
+            func: integrated function
+
+        Returns: moment approximation and error tolerance
+        """
+
+        verbose = self.params.pop('full_output', False)
+        result = quad(func, **self.params)
+        if verbose:
+            value, error, message = result
+        else:
+            value, error = result
+            message = None
+        return IntegrationResult(value, error, message)
diff --git a/src/algorithms/support_algorithms/rqmc.py b/src/algorithms/support_algorithms/rqmc.py
index 4b10c9d..543d2cb 100644
--- a/src/algorithms/support_algorithms/rqmc.py
+++ b/src/algorithms/support_algorithms/rqmc.py
@@ -5,6 +5,8 @@
 import scipy
 from numba import njit
 
+from src.algorithms.support_algorithms.integrator import IntegrationResult
+
 BITS = 30
 """Number of bits in XOR. Should be less than 64"""
 NUMBA_FAST_MATH = True
@@ -126,7 +128,8 @@ def _update(
 
         Returns:Updated mean of all rows
 
-        """
+
+                """
         values = []
         sum_of_new: float = 0.0
         for i in range(self.count):
@@ -212,9 +215,9 @@ def _xor_float(a: float, b: float) -> float:
         Returns: XOR float value
 
         """
-        a = int(a * (2**BITS))
-        b = int(b * (2**BITS))
-        return np.bitwise_xor(a, b) / 2**BITS
+        a = int(a * (2 ** BITS))
+        b = int(b * (2 ** BITS))
+        return np.bitwise_xor(a, b) / 2 ** BITS
 
     def __call__(self) -> tuple[float, float]:
         """Interface for users
@@ -223,3 +226,51 @@ def __call__(self) -> tuple[float, float]:
 
         """
         return self.rqmc()
+
+
+class RQMCIntegrator:
+    """
+    Randomize Quasi Monte Carlo Method
+
+    Args:
+        error_tolerance: pre-specified error tolerance
+        count: number of rows of random values matrix
+        base_n: number of columns of random values matrix
+        i_max: allowed number of cycles
+        a: parameter for quantile of normal distribution
+
+    """
+
+    def __init__(
+            self,
+            error_tolerance: float = 1e-6,
+            count: int = 25,
+            base_n: int = 2 ** 6,
+            i_max: int = 100,
+            a: float = 0.00047,
+    ):
+        self.error_tolerance = error_tolerance
+        self.count = count
+        self.base_n = base_n
+        self.i_max = i_max
+        self.a = a
+
+    def compute(self, func: Callable) -> IntegrationResult:
+        """
+        Compute integral via RQMC integrator
+
+        Args:
+            func: integrated function
+
+        Returns: moment approximation and error tolerance
+        """
+        result = RQMC(
+            func,
+            error_tolerance=self.error_tolerance,
+            count=self.count,
+            base_n=self.base_n,
+            i_max=self.i_max,
+            a=self.a,
+        )()
+        return IntegrationResult(result[0], result[1])
+
diff --git a/src/mixtures/abstract_mixture.py b/src/mixtures/abstract_mixture.py
index 13c8e42..3e95de4 100644
--- a/src/mixtures/abstract_mixture.py
+++ b/src/mixtures/abstract_mixture.py
@@ -5,6 +5,7 @@
 from scipy.stats import rv_continuous
 from scipy.stats.distributions import rv_frozen
 
+from src.algorithms.support_algorithms.integrator import Integrator
 
 class AbstractMixtures(metaclass=ABCMeta):
     """Base class for Mixtures"""
@@ -28,13 +29,13 @@ def __init__(self, mixture_form: str, **kwargs: Any) -> None:
             raise AssertionError(f"Unknown mixture form: {mixture_form}")
 
     @abstractmethod
-    def compute_moment(self, n: int, params: dict) -> tuple[float, float]: ...
+    def compute_moment(self, n: int, integrator: Integrator) -> tuple[float, float]: ...
 
     @abstractmethod
-    def compute_cdf(self, x: float, params: dict) -> tuple[float, float]: ...
+    def compute_cdf(self, x: float, integrator: Integrator) -> tuple[float, float]: ...
 
     @abstractmethod
-    def compute_pdf(self, x: float, params: dict) -> tuple[float, float]: ...
+    def compute_pdf(self, x: float, integrator: Integrator) -> tuple[float, float]: ...
 
     @abstractmethod
     def compute_logpdf(self, x: float, params: dict) -> tuple[float, float]: ...
diff --git a/src/mixtures/nm_mixture.py b/src/mixtures/nm_mixture.py
index 33c7580..617c498 100644
--- a/src/mixtures/nm_mixture.py
+++ b/src/mixtures/nm_mixture.py
@@ -2,13 +2,14 @@
 from typing import Any
 
 import numpy as np
-from scipy.integrate import quad
 from scipy.special import binom
 from scipy.stats import norm, rv_continuous
 from scipy.stats.distributions import rv_frozen
 
 from src.algorithms.support_algorithms.log_rqmc import LogRQMC
-from src.algorithms.support_algorithms.rqmc import RQMC
+from src.algorithms.support_algorithms.integrator import Integrator
+from src.algorithms.support_algorithms.rqmc import RQMCIntegrator
+from src.algorithms.support_algorithms.quad_integrator import QuadIntegrator
 from src.mixtures.abstract_mixture import AbstractMixtures
 
 
@@ -59,156 +60,154 @@ def _params_validation(self, data_collector: Any, params: dict[str, float | rv_c
             raise ValueError("Gamma cant be zero")
         return data_class
 
-    def _classical_moment(self, n: int, params: dict) -> tuple[float, float]:
+    def _classical_moment(self, n: int, integrator: Integrator = QuadIntegrator()) -> tuple[float, float]:
         """
         Compute n-th moment of classical NMM
 
         Args:
             n (): Moment ordinal
-            params (): Parameters of integration algorithm
+            integrator (): class of integrator with params to computing
 
         Returns: moment approximation and error tolerance
 
         """
-        mixture_moment = 0
-        error_tolerance = 0
+        mixture_moment = 0.0
+        error_tolerance = 0.0
         for k in range(0, n + 1):
             for l in range(0, k + 1):
-                coefficient = binom(n, n - k) * binom(k, k - l) * (self.params.beta ** (k - l)) * (self.params.gamma**l)
-                mixing_moment = quad(lambda u: self.params.distribution.ppf(u) ** (k - l), 0, 1, **params)
-                error_tolerance += (self.params.beta ** (k - l)) * mixing_moment[1]
-                mixture_moment += coefficient * (self.params.alpha ** (n - k)) * mixing_moment[0] * norm.moment(l)
+                coefficient = binom(n, n - k) * binom(k, k - l) * (self.params.beta ** (k - l)) * (
+                            self.params.gamma ** l)
+                mixing_moment = integrator.compute(lambda u: self.params.distribution.ppf(u) ** (k - l))
+                error_tolerance += (self.params.beta ** (k - l)) * mixing_moment.error
+                mixture_moment += coefficient * (self.params.alpha ** (n - k)) * mixing_moment.value * norm.moment(l)
         return mixture_moment, error_tolerance
 
-    def _canonical_moment(self, n: int, params: dict) -> tuple[float, float]:
+    def _canonical_moment(self, n: int, integrator: Integrator = QuadIntegrator()) -> tuple[float, float]:
         """
         Compute n-th moment of canonical NMM
 
         Args:
             n (): Moment ordinal
-            params (): Parameters of integration algorithm
+            integrator (): class of integrator with params to computing
 
         Returns: moment approximation and error tolerance
 
         """
-        mixture_moment = 0
-        error_tolerance = 0
+        mixture_moment = 0.0
+        error_tolerance = 0.0
         for k in range(0, n + 1):
-            coefficient = binom(n, n - k) * (self.params.sigma**k)
-            mixing_moment = quad(lambda u: self.params.distribution.ppf(u) ** (n - k), 0, 1, **params)
-            error_tolerance += mixing_moment[1]
-            mixture_moment += coefficient * mixing_moment[0] * norm.moment(k)
+            coefficient = binom(n, n - k) * (self.params.sigma ** k)
+            mixing_moment = integrator.compute(lambda u: self.params.distribution.ppf(u) ** (n - k))
+            error_tolerance += mixing_moment.error
+            mixture_moment += coefficient * mixing_moment.value * norm.moment(k)
         return mixture_moment, error_tolerance
 
-    def compute_moment(self, n: int, params: dict) -> tuple[float, float]:
+    def compute_moment(self, n: int, integrator: Integrator = QuadIntegrator()) -> tuple[float, float]:
         """
         Compute n-th moment of  NMM
 
         Args:
             n (): Moment ordinal
-            params (): Parameters of integration algorithm
+            integrator (): class of integrator with params to computing
 
         Returns: moment approximation and error tolerance
 
         """
         if isinstance(self.params, _NMMClassicDataCollector):
-            return self._classical_moment(n, params)
-        return self._canonical_moment(n, params)
+            return self._classical_moment(n, integrator)
+        return self._canonical_moment(n, integrator)
 
-    def _canonical_compute_cdf(self, x: float, params: dict) -> tuple[float, float]:
+    def _canonical_compute_cdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         """
         Equation for canonical cdf
         Args:
             x (): point
-            params (): parameters of RQMC algorithm
+            integrator (): class of integrator with params to computing
 
         Returns: computed cdf and error tolerance
 
         """
-        rqmc = RQMC(lambda u: norm.cdf((x - self.params.distribution.ppf(u)) / np.abs(self.params.sigma)), **params)
-        return rqmc()
+        result = integrator.compute(func=lambda u: norm.cdf((x - self.params.distribution.ppf(u)) / np.abs(self.params.sigma)))
+        return result.value, result.error
 
-    def _classical_compute_cdf(self, x: float, params: dict) -> tuple[float, float]:
+    def _classical_compute_cdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         """
         Equation for classic cdf
         Args:
             x (): point
-            params (): parameters of RQMC algorithm
+            integrator (): class of integrator with params to computing
 
         Returns: computed cdf and error tolerance
 
         """
-        rqmc = RQMC(
+        result = integrator.compute(func=
             lambda u: norm.cdf(
                 (x - self.params.alpha - self.params.beta * self.params.distribution.ppf(u)) / np.abs(self.params.gamma)
-            ),
-            **params
+            )
         )
-        return rqmc()
+        return result.value, result.error
 
-    def compute_cdf(self, x: float, params: dict) -> tuple[float, float]:
+    def compute_cdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         """
         Choose equation for cdf estimation depends on Mixture form
         Args:
             x (): point
-            params (): parameters of RQMC algorithm
+            integrator (): class of integrator with params to computing
 
         Returns: Computed pdf and error tolerance
 
         """
         if isinstance(self.params, _NMMCanonicalDataCollector):
-            return self._canonical_compute_cdf(x, params)
-        return self._classical_compute_cdf(x, params)
+            return self._canonical_compute_cdf(x, integrator)
+        return self._classical_compute_cdf(x, integrator)
 
-    def _canonical_compute_pdf(self, x: float, params: dict) -> tuple[float, float]:
+    def _canonical_compute_pdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         """
         Equation for canonical pdf
         Args:
             x (): point
-            params (): parameters of RQMC algorithm
+            integrator (): class of integrator with params to computing
 
         Returns: computed pdf and error tolerance
 
         """
-        rqmc = RQMC(
+        result = integrator.compute(func=
             lambda u: (1 / np.abs(self.params.sigma))
             * norm.pdf((x - self.params.distribution.ppf(u)) / np.abs(self.params.sigma)),
-            **params
         )
-        return rqmc()
+        return result.value, result.error
 
-    def _classical_compute_pdf(self, x: float, params: dict) -> tuple[float, float]:
+    def _classical_compute_pdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         """
         Equation for classic pdf
         Args:
             x (): point
-            params (): parameters of RQMC algorithm
+            integrator (): class of integrator with params to computing
 
         Returns: computed pdf and error tolerance
 
         """
-        rqmc = RQMC(
+        result = integrator.compute(func=
             lambda u: (1 / np.abs(self.params.gamma))
             * norm.pdf(
                 (x - self.params.alpha - self.params.beta * self.params.distribution.ppf(u)) / np.abs(self.params.gamma)
             ),
-            **params
         )
-        return rqmc()
+        return result.value, result.error
 
-    def compute_pdf(self, x: float, params: dict) -> tuple[float, float]:
+    def compute_pdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         """
         Choose equation for pdf estimation depends on Mixture form
         Args:
             x (): point
-            params (): parameters of RQMC algorithm
+            integrator (): class of integrator with params to computing
 
         Returns: Computed pdf and error tolerance
 
         """
         if isinstance(self.params, _NMMCanonicalDataCollector):
-            return self._canonical_compute_pdf(x, params)
-        return self._classical_compute_pdf(x, params)
+            return self._canonical_compute_pdf(x, integrator)
+        return self._classical_compute_pdf(x, integrator)
 
     def _classical_compute_log_pdf(self, x: float, params: dict) -> tuple[float, float]:
         """
@@ -258,4 +257,4 @@ def compute_logpdf(self, x: float, params: dict) -> tuple[float, float]:
         """
         if isinstance(self.params, _NMMCanonicalDataCollector):
             return self._canonical_compute_log_pdf(x, params)
-        return self._classical_compute_log_pdf(x, params)
+        return self._classical_compute_log_pdf(x, params)
\ No newline at end of file
diff --git a/src/mixtures/nmv_mixture.py b/src/mixtures/nmv_mixture.py
index c090ff3..a147cba 100644
--- a/src/mixtures/nmv_mixture.py
+++ b/src/mixtures/nmv_mixture.py
@@ -8,7 +8,8 @@
 from scipy.stats.distributions import rv_frozen
 
 from src.algorithms.support_algorithms.log_rqmc import LogRQMC
-from src.algorithms.support_algorithms.rqmc import RQMC
+from src.algorithms.support_algorithms.integrator import Integrator
+from src.algorithms.support_algorithms.rqmc import RQMCIntegrator
 from src.mixtures.abstract_mixture import AbstractMixtures
 
 
@@ -40,13 +41,13 @@ class NormalMeanVarianceMixtures(AbstractMixtures):
     def __init__(self, mixture_form: str, **kwargs: Any) -> None:
         super().__init__(mixture_form, **kwargs)
 
-    def _classical_moment(self, n: int, params: dict) -> tuple[float, float]:
+    def _classical_moment(self, n: int, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         """
         Compute n-th moment of classical NMM
 
         Args:
             n (): Moment ordinal
-            params (): Parameters of integration algorithm
+            integrator (): class of integrator with params to computing
 
         Returns: moment approximation and error tolerance
 
@@ -67,16 +68,16 @@ def integral_func(u: float) -> float:
                     )
             return result
 
-        rqmc = RQMC(lambda u: integral_func(u), **params)
-        return rqmc()
+        result = integrator.compute(func=lambda u: integral_func(u))
+        return result.value, result.error
 
-    def _canonical_moment(self, n: int, params: dict) -> tuple[float, float]:
+    def _canonical_moment(self, n: int, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         """
         Compute n-th moment of classical NMM
 
         Args:
             n (): Moment ordinal
-            params (): Parameters of integration algorithm
+            integrator (): class of integrator with params to computing
 
         Returns: moment approximation and error tolerance
 
@@ -96,15 +97,15 @@ def integral_func(u: float) -> float:
                     )
             return result
 
-        rqmc = RQMC(lambda u: integral_func(u), **params)
-        return rqmc()
+        result = integrator.compute(func=lambda u: integral_func(u))
+        return result.value, result.error
 
-    def compute_moment(self, n: int, params: dict) -> tuple[float, float]:
+    def compute_moment(self, n: int, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         if isinstance(self.params, _NMVMClassicDataCollector):
-            return self._classical_moment(n, params)
-        return self._canonical_moment(n, params)
+            return self._classical_moment(n, integrator)
+        return self._canonical_moment(n, integrator)
 
-    def _classical_cdf(self, x: float, params: dict) -> tuple[float, float]:
+    def _classical_cdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         def _inner_func(u: float) -> float:
             ppf = lru_cache()(self.params.distribution.ppf)(u)
             point = (x - self.params.alpha) / (np.sqrt(ppf) * self.params.gamma) - (
@@ -112,24 +113,24 @@ def _inner_func(u: float) -> float:
             )
             return norm.cdf(point)
 
-        rqmc = RQMC(lambda u: _inner_func(u), **params)
-        return rqmc()
+        result = integrator.compute(func=lambda u: _inner_func(u))
+        return result.value, result.error
 
-    def _canonical_cdf(self, x: float, params: dict) -> tuple[float, float]:
+    def _canonical_cdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         def _inner_func(u: float) -> float:
             ppf = self.params.distribution.ppf(u)
             point = (x - self.params.alpha) / (np.sqrt(ppf)) - (self.params.mu * np.sqrt(ppf))
             return norm.cdf(point)
 
-        rqmc = RQMC(lambda u: _inner_func(u), **params)
-        return rqmc()
+        result = integrator.compute(func=lambda u: _inner_func(u))
+        return result.value, result.error
 
-    def compute_cdf(self, x: float, params: dict) -> tuple[float, float]:
+    def compute_cdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         if isinstance(self.params, _NMVMClassicDataCollector):
-            return self._classical_cdf(x, params)
-        return self._canonical_cdf(x, params)
+            return self._classical_cdf(x, integrator)
+        return self._canonical_cdf(x, integrator)
 
-    def _classical_pdf(self, x: float, params: dict) -> tuple[float, float]:
+    def _classical_pdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         def _inner_func(u: float) -> float:
             ppf = self.params.distribution.ppf(u)
             return (
@@ -140,10 +141,10 @@ def _inner_func(u: float) -> float:
                 )
             )
 
-        rqmc = RQMC(lambda u: _inner_func(u), **params)()
-        return np.exp(self.params.beta * (x - self.params.alpha) / self.params.gamma**2) * rqmc[0], rqmc[1]
+        result = integrator.compute(func=lambda u: _inner_func(u))
+        return np.exp(self.params.beta * (x - self.params.alpha) / self.params.gamma**2) * result.value, result.error
 
-    def _canonical_pdf(self, x: float, params: dict) -> tuple[float, float]:
+    def _canonical_pdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         def _inner_func(u: float) -> float:
             ppf = self.params.distribution.ppf(u)
             return (
@@ -152,9 +153,8 @@ def _inner_func(u: float) -> float:
                 * np.exp(-((x - self.params.alpha) ** 2 + self.params.mu**2 * ppf**2) / (2 * ppf))
             )
 
-        rqmc = RQMC(lambda u: _inner_func(u), **params)
-        res = rqmc()
-        return np.exp(self.params.mu * (x - self.params.alpha)) * res[0], res[1]
+        result = integrator.compute(func=lambda u: _inner_func(u))
+        return np.exp(self.params.mu * (x - self.params.alpha)) * result.value, result.error
 
     def _classical_log_pdf(self, x: float, params: dict) -> tuple[float, float]:
         def _inner_func(u: float) -> float:
@@ -165,8 +165,8 @@ def _inner_func(u: float) -> float:
                 + ppf * self.params.gamma**2 * np.log(2 * np.pi * ppf * self.params.gamma**2)
             ) / (2 * ppf * self.params.gamma**2)
 
-        rqmc = LogRQMC(lambda u: _inner_func(u), **params)
-        return rqmc()
+        result = LogRQMC(lambda u: _inner_func(u), **params)
+        return result()
 
     def _canonical_log_pdf(self, x: float, params: dict) -> tuple[float, float]:
         def _inner_func(u: float) -> float:
@@ -175,13 +175,13 @@ def _inner_func(u: float) -> float:
                 2 * ppf
             )
 
-        rqmc = LogRQMC(lambda u: _inner_func(u), **params)
-        return rqmc()
+        result = LogRQMC(lambda u: _inner_func(u), **params)
+        return result()
 
-    def compute_pdf(self, x: float, params: dict) -> tuple[float, float]:
+    def compute_pdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         if isinstance(self.params, _NMVMClassicDataCollector):
-            return self._classical_pdf(x, params)
-        return self._canonical_pdf(x, params)
+            return self._classical_pdf(x, integrator)
+        return self._canonical_pdf(x, integrator)
 
     def compute_logpdf(self, x: float, params: dict) -> tuple[Any, float]:
         if isinstance(self.params, _NMVMClassicDataCollector):
diff --git a/src/mixtures/nv_mixture.py b/src/mixtures/nv_mixture.py
index 7905e6e..ffc4f08 100644
--- a/src/mixtures/nv_mixture.py
+++ b/src/mixtures/nv_mixture.py
@@ -8,7 +8,8 @@
 from scipy.stats.distributions import rv_frozen
 
 from src.algorithms.support_algorithms.log_rqmc import LogRQMC
-from src.algorithms.support_algorithms.rqmc import RQMC
+from src.algorithms.support_algorithms.integrator import Integrator
+from src.algorithms.support_algorithms.rqmc import RQMCIntegrator
 from src.mixtures.abstract_mixture import AbstractMixtures
 
 
@@ -39,12 +40,12 @@ class NormalVarianceMixtures(AbstractMixtures):
     def __init__(self, mixture_form: str, **kwargs: Any) -> None:
         super().__init__(mixture_form, **kwargs)
 
-    def compute_moment(self, n: int, params: dict) -> tuple[float, float]:
+    def compute_moment(self, n: int, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         """
         Compute n-th moment of  NVM
         Args:
             n (): Moment ordinal
-            params (): Parameters of integration algorithm
+            integrator (): class of integrator with params to computing
         Returns: moment approximation and error tolerance
         """
         gamma = self.params.gamma if isinstance(self.params, _NVMClassicDataCollector) else 1
@@ -60,16 +61,15 @@ def integrate_func(u: float) -> float:
                     for k in range(0, n + 1)
                 ]
             )
+        result = integrator.compute(func=integrate_func)
+        return result.value, result.error
 
-        result = RQMC(integrate_func, **params)()
-        return result
-
-    def compute_cdf(self, x: float, params: dict) -> tuple[float, float]:
+    def compute_cdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         parametric_norm = norm(0, self.params.gamma if isinstance(self.params, _NVMClassicDataCollector) else 1)
-        rqmc = RQMC(
-            lambda u: parametric_norm.cdf((x - self.params.alpha) / np.sqrt(self.params.distribution.ppf(u))), **params
+        result = integrator.compute(func=
+            lambda u: parametric_norm.cdf((x - self.params.alpha) / np.sqrt(self.params.distribution.ppf(u)))
         )
-        return rqmc()
+        return result.value, result.error
 
     @lru_cache()
     def _integrand_func(self, u: float, d: float, gamma: float) -> float:
@@ -80,11 +80,11 @@ def _log_integrand_func(self, u: float, d: float, gamma: float | int | np.int64)
         ppf = self.params.distribution.ppf(u)
         return -(ppf * np.log(np.pi * 2 * ppf * gamma**2) + d) / (2 * ppf)
 
-    def compute_pdf(self, x: float, params: dict) -> tuple[float, float]:
+    def compute_pdf(self, x: float, integrator: Integrator = RQMCIntegrator()) -> tuple[float, float]:
         gamma = self.params.gamma if isinstance(self.params, _NVMClassicDataCollector) else 1
         d = (x - self.params.alpha) ** 2 / gamma**2
-        rqmc = RQMC(lambda u: self._integrand_func(u, d, gamma), **params)
-        return rqmc()
+        result = integrator.compute(func=lambda u: self._integrand_func(u, d, gamma))
+        return result.value, result.error
 
     def compute_logpdf(self, x: float, params: dict) -> tuple[float, float]:
         gamma = self.params.gamma if isinstance(self.params, _NVMClassicDataCollector) else 1