Moments method

examples/big_mono_tests.py

@@ -3,17 +3,17 @@
 import random
 
 import numpy as np
-
-from examples.config import MAX_WORKERS, TESTS_OPTIMIZERS
-from examples.mono_test_generator import generate_mono_test
-from examples.utils import Clicker, Test, init_solver, run_tests, save_results
 from mpest.models import (
     AModelWithGenerator,
     ExponentialModel,
     GaussianModel,
     WeibullModelExp,
 )
 
+from examples.config import MAX_WORKERS, TESTS_OPTIMIZERS
+from examples.mono_test_generator import generate_mono_test
+from examples.utils import Clicker, Test, init_solver, run_tests, save_results
+
 if __name__ == "__main__":
     random.seed(42)
     np.random.seed(42)

examples/diff_test.py

@@ -3,12 +3,12 @@
 import random
 
 import numpy as np
+from mpest import Distribution, MixtureDistribution, Problem
+from mpest.models import GaussianModel, WeibullModelExp
 
 from examples.config import MAX_WORKERS, TESTS_OPTIMIZERS
 from examples.mono_test_generator import Clicker
 from examples.utils import Test, init_solver, run_tests, save_results
-from mpest import Distribution, MixtureDistribution, Problem
-from mpest.models import GaussianModel, WeibullModelExp
 
 # Gaussian
 

examples/mono_test_generator.py

@@ -4,14 +4,14 @@
 from collections.abc import Iterable
 
 import numpy as np
-
-from examples.utils import Clicker, Test
 from mpest.core.distribution import Distribution
 from mpest.core.mixture_distribution import MixtureDistribution
 from mpest.core.problem import Problem
 from mpest.em import EM
 from mpest.models import AModel, AModelWithGenerator
 
+from examples.utils import Clicker, Test
+
 
 def generate_mono_test(
     model_t: type[AModelWithGenerator],

examples/prepare_result.py

@@ -4,13 +4,13 @@
 
 import numpy as np
 import pandas as pd
+from mpest.annotations import Samples
+from mpest.core.mixture_distribution import DistributionInMixture, MixtureDistribution
 from tqdm.contrib.concurrent import process_map
 
 from examples.config import MAX_WORKERS
 from examples.mono_test_generator import Clicker
 from examples.utils import SingleSolverResult, TestResult
-from mpest.annotations import Samples
-from mpest.core.mixture_distribution import DistributionInMixture, MixtureDistribution
 
 
 def nll(samples: Samples, mixture: MixtureDistribution) -> float:

examples/quick_test.py

@@ -3,10 +3,6 @@
 import random
 
 import numpy as np
-
-from examples.config import MAX_WORKERS
-from examples.mono_test_generator import generate_mono_test
-from examples.utils import Clicker, Test, init_solver, run_tests, save_results
 from mpest.models import (
     AModelWithGenerator,
     ExponentialModel,
@@ -15,6 +11,10 @@
 )
 from mpest.optimizers import ALL_OPTIMIZERS
 
+from examples.config import MAX_WORKERS
+from examples.mono_test_generator import generate_mono_test
+from examples.utils import Clicker, Test, init_solver, run_tests, save_results
+
 
 def run_test():
     """Runs the mixture distributions of single model quick test"""

examples/readme_example/example.py

@@ -3,7 +3,6 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import seaborn as sns
-
 from mpest import Distribution, MixtureDistribution, Problem
 from mpest.em import EM
 from mpest.em.breakpointers import StepCountBreakpointer

examples/readme_example/example_ml.py

@@ -7,11 +7,6 @@
 import numpy as np
 import pandas as pd
 import seaborn as sns
-from scipy.stats import entropy, wasserstein_distance
-from sklearn.cluster import DBSCAN, AgglomerativeClustering, KMeans
-from sklearn.metrics import calinski_harabasz_score, davies_bouldin_score, silhouette_score
-from sklearn.neighbors import NearestNeighbors
-
 from mpest import Distribution, MixtureDistribution, Problem
 from mpest.em import EM
 from mpest.em.breakpointers import StepCountBreakpointer
@@ -20,6 +15,10 @@
 from mpest.em.methods.method import Method
 from mpest.models import GaussianModel, WeibullModelExp
 from mpest.optimizers import ScipyCG
+from scipy.stats import entropy, wasserstein_distance
+from sklearn.cluster import DBSCAN, AgglomerativeClustering, KMeans
+from sklearn.metrics import calinski_harabasz_score, davies_bouldin_score, silhouette_score
+from sklearn.neighbors import NearestNeighbors
 
 os.makedirs("results", exist_ok=True)
 os.makedirs("results/plots", exist_ok=True)

examples/utils.py

@@ -7,9 +7,6 @@
 from typing import ClassVar, NamedTuple
 
 import numpy as np
-from tqdm.contrib.concurrent import process_map
-
-from examples.config import RESULTS_FOLDER
 from mpest.annotations import Samples
 from mpest.core.mixture_distribution import MixtureDistribution
 from mpest.core.problem import Problem, Result
@@ -21,6 +18,9 @@
 )
 from mpest.em.methods.likelihood_method import LikelihoodMethod
 from mpest.optimizers import TOptimizer
+from tqdm.contrib.concurrent import process_map
+
+from examples.config import RESULTS_FOLDER
 
 np.seterr(all="ignore")
 

mpest/em/methods/moments_method.py

@@ -0,0 +1,82 @@
+"""The module in which the moments method is presented"""
+
+import numpy as np
+
+from mpest import Samples
+from mpest.core.distribution import Distribution
+from mpest.core.mixture_distribution import MixtureDistribution
+from mpest.core.problem import Problem, Result
+from mpest.em.methods.abstract_steps import AMaximization
+from mpest.exceptions import MStepError
+from mpest.utils import ResultWithError
+
+EResult = tuple[Problem, np.ndarray] | ResultWithError[MixtureDistribution]
+
+
+class MomentsMStep(AMaximization[EResult]):
+    """
+    Class which calculate new params using matrix with indicator from E step.
+    """
+
+    def calc_order_moment_of_index_element(self, order: int, i: int, samples: Samples, indicators: np.ndarray) -> float:
+        """
+        A function that calculates the list of n-th moments of each distribution.
+
+        :param order: Order of Moment.
+        :param i: The number of the distribution for which we count the moment.
+        :param samples: Ndarray with samples.
+        :param indicators: Matrix with indicators
+
+        :return:  order-Moment of index element.
+        """
+
+        sum_j_row_probabilities = np.sum(indicators[i])
+
+        if sum_j_row_probabilities == 0:
+            return 0
+
+        moment_values = samples**order
+
+        numerator = np.sum(moment_values * indicators[i])
+
+        return numerator / sum_j_row_probabilities
+
+    def step(self, e_result: EResult) -> Result:
+        """
+        A function that performs M step
+
+        :param e_result: Tuple with problem, new_priors and indicators.
+        """
+
+        if isinstance(e_result, ResultWithError):
+            return e_result
+
+        problem, indicators = e_result
+
+        samples = problem.samples
+
+        mixture = problem.distributions
+
+        new_priors = np.sum(indicators, axis=1) / len(samples)
+
+        max_params_count = max(len(d.params) for d in mixture)
+        moments = np.zeros(shape=[len(mixture), max_params_count])
+
+        for j, d in enumerate(mixture):
+            for r in range(len(d.params)):
+                moments[j][r] = self.calc_order_moment_of_index_element(r + 1, j, samples, indicators)
+
+        for i, d in enumerate(mixture):
+            if d.model.name == "WeibullExp" and (moments[i][0] * moments[i][1] < 0):
+                error = MStepError("The Weibull distribution degenerated in the first step.")
+                return ResultWithError(mixture.distributions, error)
+
+        new_distributions = []
+
+        for j, d in enumerate(mixture):
+            new_params = d.model.calc_moments_params(moments[j])
+            new_d = Distribution(d.model, d.model.params_convert_to_model(new_params))
+            new_distributions.append(new_d)
+
+        new_mixture = MixtureDistribution.from_distributions(new_distributions, new_priors)
+        return ResultWithError(new_mixture)

mpest/models/exponential.py

@@ -65,3 +65,13 @@ def calc_params(self, moments: list[float]):
         lm = 1 / moments[0]
 
         return np.array([lm])
+
+    def calc_moments_params(self, moments: list[float]):
+        """
+        The function for calculating params using moments
+        """
+
+        # Calculate lambda parameter
+        lm = 1 / moments[0]
+
+        return np.array([lm])

mpest/models/gaussian.py

@@ -74,3 +74,21 @@ def calc_params(self, moments: list[float]) -> np.ndarray:
         variance = m2 * np.sqrt(np.pi)
 
         return np.array([mean, variance])
+
+    def calc_moments_params(self, moments: list[float]) -> np.ndarray:
+        """
+        The function for calculating params using moments
+        """
+
+        m1 = moments[0]
+        m2 = moments[1]
+
+        # Calculate mean parameter
+        mu = m1
+
+        # Calculate variance parameter
+        variance = m2 - m1**2
+
+        sigma = np.sqrt(variance)
+
+        return np.array([mu, sigma])

mpest/models/weibull.py

@@ -3,6 +3,8 @@
 import math
 
 import numpy as np
+from scipy.optimize import root_scalar
+from scipy.special import gamma
 from scipy.stats import weibull_min
 
 from mpest.annotations import Params, Samples
@@ -86,3 +88,25 @@ def calc_params(self, moments: list[float]):
         lm = m1 / math.gamma(1 + 1 / k)
 
         return np.array([k, lm])
+
+    def calc_moments_params(self, moments: list[float]):
+        """
+        The function for calculating params using moments
+        """
+
+        m1, m2 = moments[0], moments[1]
+
+        moments_ratio = m2 / (m1**2)
+
+        def equation_for_k(k):
+            return gamma(1 + 2 / k) / (gamma(1 + 1 / k) ** 2) - moments_ratio
+
+        solution = root_scalar(equation_for_k, method="brentq", bracket=[0.02, 100])
+        if not solution.converged:
+            raise RuntimeError(f"Error in calculating the equation: m1={m1}, m2={m2}")
+
+        k = solution.root
+
+        lm = m1 / gamma(1 + 1 / k)
+
+        return np.array([k, lm])

tests/tests_moments/moments_utils.py

@@ -0,0 +1,21 @@
+from mpest.core.problem import Problem, Result
+from mpest.em import EM
+from mpest.em.breakpointers import ParamDifferBreakpointer, StepCountBreakpointer
+from mpest.em.distribution_checkers import (
+    FiniteChecker,
+    PriorProbabilityThresholdChecker,
+)
+from mpest.em.methods.likelihood_method import BayesEStep
+from mpest.em.methods.method import Method
+from mpest.em.methods.moments_method import MomentsMStep
+
+
+def run_test(problem: Problem, deviation: float) -> Result:
+    method = Method(BayesEStep(), MomentsMStep())
+    em_algo = EM(
+        StepCountBreakpointer() + ParamDifferBreakpointer(deviation=deviation),
+        FiniteChecker() + PriorProbabilityThresholdChecker(),
+        method,
+    )
+
+    return em_algo.solve(problem=problem)