diff --git a/divik/_matlab_legacy.py b/divik/_matlab_legacy.py
index b08f1c72..7b793562 100644
--- a/divik/_matlab_legacy.py
+++ b/divik/_matlab_legacy.py
@@ -45,3 +45,35 @@ def find_thresholds(values: np.ndarray, max_components: int = 10) -> np.ndarray:
         return np.array([])
     with sigsegv_handler():
         return gn.find_thresholds(values - offset, max_components) + offset
+
+def find_gaussian_mixtures(x: np.ndarray, counts:np.ndarray, KS: int = 10) -> dict:
+    """Find parameters for decomposition of values by GMM.
+
+    Parameters
+    ----------
+    values:
+        vector of values to decompose
+
+    KS:
+        number of Gaussian components
+
+    Returns
+    -------
+    dictionary with arrays of weights, means, standard deviations and log likelihood of decomposition
+    """
+    if x.size <= KS:
+        KS = x.size
+    if x.size == 0:
+        return {"weights": np.array([]), "mu": np.array([]), "sig": np.array([]), "TIC": np.nan, "l_lik": np.nan}
+    if KS <= 0:
+        raise ValueError("KS must be positive")
+    if x.ndim != 1:
+        raise ValueError("x must be 1D array")
+    if counts.ndim != 1:
+        raise ValueError("counts must be 1D array")
+    x = np.ascontiguousarray(x)
+    counts = np.ascontiguousarray(counts)
+    if np.min(x) == np.max(x):
+        return {"weights": np.array([]), "mu": np.array([]), "sig": np.array([]), "TIC": np.nan, "l_lik": np.nan}
+    with sigsegv_handler():
+        return gn.find_gaussian_mixtures(x, counts, KS)
diff --git a/gamred_native/fetch_thresholds_data.c b/gamred_native/fetch_thresholds_data.c
index e4d6bb33..e864c5f6 100644
--- a/gamred_native/fetch_thresholds_data.c
+++ b/gamred_native/fetch_thresholds_data.c
@@ -14,5 +14,5 @@
 /* Variable Definitions */
 omp_nest_lock_t emlrtNestLockGlobal;
 boolean_T isInitialized_fetch_thresholds = false;
-
+boolean_T isInitialized_gaussian_mixture_simple = false;
 /* End of code generation (fetch_thresholds_data.c) */
diff --git a/gamred_native/fetch_thresholds_data.h b/gamred_native/fetch_thresholds_data.h
index f712e83d..d3e459a3 100644
--- a/gamred_native/fetch_thresholds_data.h
+++ b/gamred_native/fetch_thresholds_data.h
@@ -19,7 +19,7 @@
 /* Variable Declarations */
 extern omp_nest_lock_t emlrtNestLockGlobal;
 extern boolean_T isInitialized_fetch_thresholds;
-
+extern boolean_T isInitialized_gaussian_mixture_simple;
 #endif
 
 /* End of code generation (fetch_thresholds_data.h) */
diff --git a/gamred_native/gaussian_mixture_simple.c b/gamred_native/gaussian_mixture_simple.c
index 5cbfd544..a9c8ffac 100644
--- a/gamred_native/gaussian_mixture_simple.c
+++ b/gamred_native/gaussian_mixture_simple.c
@@ -7,10 +7,15 @@
  */
 
 /* Include files */
+#include <stdio.h>
+#include <stdlib.h>
+
 #include "gaussian_mixture_simple.h"
 #include "dyn_pr_split_w.h"
 #include "fetch_thresholds.h"
+#include "fetch_thresholds_data.h"
 #include "fetch_thresholds_emxutil.h"
+#include "gaussian_mixture_simple_initialize.h"
 #include "g_mix_est_fast_lik.h"
 #include "histcounts.h"
 #include "mean.h"
@@ -57,7 +62,11 @@ void gaussian_mixture_simple(const emxArray_real_T *x, const emxArray_real_T
   double d;
   double wwec_tmp;
   bool exitg1;
+
   emxInit_real_T(&b_x, 1);
+  if (isInitialized_gaussian_mixture_simple == false) {
+    gaussian_mixture_simple_initialize();
+  }
 
   /*  FUNCTION:  */
   /*  gaussian_mixture */
@@ -80,6 +89,7 @@ void gaussian_mixture_simple(const emxArray_real_T *x, const emxArray_real_T
   /*  Andrzej Polanski */
   /*  email: andrzej.polanski@polsl.pl */
   /* 'gaussian_mixture_simple:29' TIC=sum(x.*counts); */
+  
   i = b_x->size[0];
   b_x->size[0] = x->size[0];
   emxEnsureCapacity_real_T(b_x, i);
@@ -99,9 +109,9 @@ void gaussian_mixture_simple(const emxArray_real_T *x, const emxArray_real_T
   }
 
   emxFree_real_T(&b_x);
-
+  
   /* 'gaussian_mixture_simple:30' if KS==1 */
-  if (KS == 1.0) {
+  if (KS == 1) {
     /* 'gaussian_mixture_simple:31' mu_est=nanmean(x); */
     if (x->size[0] == 0) {
       y = rtNaN;
@@ -122,7 +132,7 @@ void gaussian_mixture_simple(const emxArray_real_T *x, const emxArray_real_T
         y /= (double)idx;
       }
     }
-
+    
     i = mu_est->size[0];
     mu_est->size[0] = 1;
     emxEnsureCapacity_real_T(mu_est, i);
@@ -143,6 +153,7 @@ void gaussian_mixture_simple(const emxArray_real_T *x, const emxArray_real_T
     /* 'gaussian_mixture_simple:34' l_lik=NaN; */
     *l_lik = rtNaN;
   } else {
+
     emxInit_real_T(&y_out, 2);
     emxInit_real_T(&edges, 2);
 
@@ -167,7 +178,7 @@ void gaussian_mixture_simple(const emxArray_real_T *x, const emxArray_real_T
       i = 1;
       i1 = edges->size[1];
     }
-
+    
     emxInit_real_T(&b_edges, 2);
     i2 = b_edges->size[0] * b_edges->size[1];
     b_edges->size[0] = 2;
@@ -176,12 +187,12 @@ void gaussian_mixture_simple(const emxArray_real_T *x, const emxArray_real_T
     for (i2 = 0; i2 < loop_ub; i2++) {
       b_edges->data[2 * i2] = edges->data[i2];
     }
-
+    
     loop_ub = i1 - i;
     for (i1 = 0; i1 < loop_ub; i1++) {
       b_edges->data[2 * i1 + 1] = edges->data[i + i1];
     }
-
+    
     emxFree_real_T(&edges);
     emxInit_real_T(&mz_out, 2);
     emxInit_real_T(&aux_mx, 2);
@@ -196,6 +207,7 @@ void gaussian_mixture_simple(const emxArray_real_T *x, const emxArray_real_T
     /* 'dyn_pr_split_w_aux:8' N=length(data); */
     /*  aux_mx */
     /* 'dyn_pr_split_w_aux:11' aux_mx=zeros(N,N); */
+    
     i = aux_mx->size[0] * aux_mx->size[1];
     aux_mx->size[0] = mz_out->size[1];
     aux_mx->size[1] = mz_out->size[1];
@@ -205,11 +217,12 @@ void gaussian_mixture_simple(const emxArray_real_T *x, const emxArray_real_T
     for (i = 0; i < loop_ub; i++) {
       aux_mx->data[i] = 0.0;
     }
-
+    
     /* 'dyn_pr_split_w_aux:12' for kk=1:N-1 */
     i = mz_out->size[1];
     emxInit_real_T(&b_mz_out, 2);
     emxInit_real_T(&b_y_out, 2);
+    
     for (idx = 0; idx <= i - 2; idx++) {
       /* 'dyn_pr_split_w_aux:13' for jj=kk+1:N */
       i1 = mz_out->size[1] - idx;
@@ -259,6 +272,7 @@ void gaussian_mixture_simple(const emxArray_real_T *x, const emxArray_real_T
     emxInit_real_T(&Q, 2);
 
     /* 'gaussian_mixture_simple:42' [Q,opt_part]=dyn_pr_split_w(mz_out,y_out,KS-1,aux_mx); */
+
     dyn_pr_split_w(mz_out, y_out, KS - 1.0, aux_mx, Q, pp_ini);
 
     /* 'gaussian_mixture_simple:43' part_cl=[1 opt_part Nb+1]; */
@@ -489,6 +503,7 @@ void gaussian_mixture_simple(const emxArray_real_T *x, const emxArray_real_T
     emxFree_real_T(&mu_ini);
     emxFree_real_T(&pp_ini);
   }
+
 }
 
 /* End of code generation (gaussian_mixture_simple.c) */
diff --git a/gamred_native/gaussian_mixture_simple_initialize.c b/gamred_native/gaussian_mixture_simple_initialize.c
new file mode 100644
index 00000000..bb4c8a08
--- /dev/null
+++ b/gamred_native/gaussian_mixture_simple_initialize.c
@@ -0,0 +1,23 @@
+/*
+ *
+ * gaussian_mixture_simple_initialize.c
+ *
+ * Code generation for function 'gaussian_mixture_simple_initialize'
+ *
+ */
+
+/* Include files */
+#include "gaussian_mixture_simple_initialize.h"
+#include "fetch_thresholds.h"
+#include "fetch_thresholds_data.h"
+#include "rt_nonfinite.h"
+
+/* Function Definitions */
+void gaussian_mixture_simple_initialize(void)
+{
+  rt_InitInfAndNaN();
+  omp_init_nest_lock(&emlrtNestLockGlobal);
+  isInitialized_gaussian_mixture_simple = true;
+}
+
+/* End of code generation (gaussian_mixture_simple_initialize.c) */
diff --git a/gamred_native/gaussian_mixture_simple_initialize.h b/gamred_native/gaussian_mixture_simple_initialize.h
new file mode 100644
index 00000000..0ea750d9
--- /dev/null
+++ b/gamred_native/gaussian_mixture_simple_initialize.h
@@ -0,0 +1,24 @@
+/*
+ *
+ * gaussian_mixture_simple_initialize.h
+ *
+ * Code generation for function 'gaussian_mixture_simple_initialize'
+ *
+ */
+
+#ifndef GAUSSIAN_MIXTURE_SIMPLE_INITIALIZE_H
+#define GAUSSIAN_MIXTURE_SIMPLE_INITIALIZE_H
+
+/* Include files */
+#include <stddef.h>
+#include <stdlib.h>
+#include "rtwtypes.h"
+#include "omp.h"
+#include "fetch_thresholds_types.h"
+
+/* Function Declarations */
+extern void gaussian_mixture_simple_initialize(void);
+
+#endif
+
+/* End of code generation (gaussian_mixture_simple_initialize.h) */
diff --git a/gamred_native/gaussian_mixture_simple_terminate.c b/gamred_native/gaussian_mixture_simple_terminate.c
new file mode 100644
index 00000000..936c3e55
--- /dev/null
+++ b/gamred_native/gaussian_mixture_simple_terminate.c
@@ -0,0 +1,22 @@
+/*
+ *
+ * fetch_thresholds_terminate.c
+ *
+ * Code generation for function 'gaussian_mixture_simple_terminate'
+ *
+ */
+
+/* Include files */
+#include "gaussian_mixture_simple_terminate.h"
+#include "fetch_thresholds.h"
+#include "fetch_thresholds_data.h"
+#include "rt_nonfinite.h"
+
+/* Function Definitions */
+void gaussian_mixture_simple_terminate(void)
+{
+  omp_destroy_nest_lock(&emlrtNestLockGlobal);
+  isInitialized_gaussian_mixture_simple = false;
+}
+
+/* End of code generation (gaussian_mixture_simple_terminate.c) */
diff --git a/gamred_native/gaussian_mixture_simple_terminate.h b/gamred_native/gaussian_mixture_simple_terminate.h
new file mode 100644
index 00000000..2212092c
--- /dev/null
+++ b/gamred_native/gaussian_mixture_simple_terminate.h
@@ -0,0 +1,24 @@
+/*
+ *
+ * gaussian_mixture_simple_terminate.h
+ *
+ * Code generation for function 'gaussian_mixture_simple_terminate'
+ *
+ */
+
+#ifndef GAUSSIAN_MIXTURE_SIMPLE_TERMINATE_H
+#define GAUSSIAN_MIXTURE_SIMPLE_TERMINATE_H
+
+/* Include files */
+#include <stddef.h>
+#include <stdlib.h>
+#include "rtwtypes.h"
+#include "omp.h"
+#include "fetch_thresholds_types.h"
+
+/* Function Declarations */
+extern void gaussian_mixture_simple_terminate(void);
+
+#endif
+
+/* End of code generation (gaussian_mixture_simple_terminate.h) */
diff --git a/gamred_native/main.c b/gamred_native/main.c
index 06c80c76..cfdb614d 100644
--- a/gamred_native/main.c
+++ b/gamred_native/main.c
@@ -37,8 +37,10 @@
 #include "fetch_thresholds.h"
 #include "main.h"
 #include "fetch_thresholds_terminate.h"
+#include "gaussian_mixture_simple_terminate.h"
 #include "fetch_thresholds_emxAPI.h"
 #include "fetch_thresholds_initialize.h"
+#include "gaussian_mixture_simple_initialize.h"
 
 static emxArray_real_T *PyArrayObject_to_emxArray(PyArrayObject *vals)
 {
@@ -56,6 +58,75 @@ static PyArrayObject *emxArray_to_PyArrayObject(emxArray_real_T *vals)
   return result;
 }
 
+static PyObject *method_find_gaussian_mixtures(PyObject *self, PyObject *args) {
+
+    PyObject *x=NULL;
+    PyArrayObject *x_arr=NULL;
+
+    PyObject *counts=NULL;
+    PyArrayObject *counts_arr=NULL;
+    
+    PyArrayObject *pp=NULL;
+    PyArrayObject *mu=NULL;
+    PyArrayObject *sig=NULL;
+    float l_lik = 0.0;
+    float TIC = 0.0;
+    unsigned long KS = 0;
+    emxArray_real_T *native_x;
+    emxArray_real_T *native_counts;
+    emxArray_real_T *native_mu;
+    emxArray_real_T *native_pp;
+    emxArray_real_T *native_sig;
+
+    /* Parse arguments */
+    // https://docs.python.org/3/c-api/arg.html
+    if(!PyArg_ParseTuple(args, "OOk", &x, &counts, &KS)) {
+
+        return NULL;
+    }
+
+    x_arr = (PyArrayObject *)PyArray_FROM_OTF(x, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
+    if (x_arr == NULL) {
+        return NULL;
+    }
+
+    counts_arr = (PyArrayObject *)PyArray_FROM_OTF(counts, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);
+    if (counts_arr == NULL) {
+        return NULL;
+    }
+
+
+    /* Initialize the application.
+      You do not need to do this more than one time. */
+    gaussian_mixture_simple_initialize();
+
+    emxInitArray_real_T(&native_pp, 1);
+    emxInitArray_real_T(&native_mu, 1);
+    emxInitArray_real_T(&native_sig, 1);
+
+    native_counts = PyArrayObject_to_emxArray(counts_arr);
+    native_x = PyArrayObject_to_emxArray(x_arr);
+    Py_DECREF(counts_arr);
+    Py_DECREF(x_arr);
+
+    gaussian_mixture_simple(native_x, native_counts, (double)KS, native_pp, native_mu, native_sig, &TIC, &l_lik);
+
+    emxDestroyArray_real_T(native_x);
+    emxDestroyArray_real_T(native_counts);
+
+    mu = emxArray_to_PyArrayObject(native_mu);
+    emxDestroyArray_real_T(native_mu);
+    sig = emxArray_to_PyArrayObject(native_sig);
+    emxDestroyArray_real_T(native_sig);
+    pp = emxArray_to_PyArrayObject(native_pp);
+    emxDestroyArray_real_T(native_pp);
+    
+    /* Terminate the application.
+      You do not need to do this more than one time. */
+    gaussian_mixture_simple_terminate();
+    return Py_BuildValue("{s:O,s:O,s:O,s:f,s:f}", "weights", (PyObject*)pp, "mu", (PyObject*)mu, "sigma", (PyObject*)sig, "TIC", TIC, "l_lik", l_lik);
+}
+
 static PyObject *method_find_thresholds(PyObject *self, PyObject *args) {
     PyObject *vals=NULL;
     PyArrayObject *vals_arr=NULL;
@@ -98,14 +169,13 @@ static PyObject *method_find_thresholds(PyObject *self, PyObject *args) {
 
 static PyMethodDef GamredNativeMethods[] = {
     {"find_thresholds", method_find_thresholds, METH_VARARGS, "Python interface for the MATLAB fetch_thresholds function"},
-    {NULL, NULL, 0, NULL}
+    {"find_gaussian_mixtures", method_find_gaussian_mixtures, METH_VARARGS, "Python interface for the MATLAB find_gaussian_mixtures function"}
 };
 
-
 static struct PyModuleDef gamred_native_module = {
     PyModuleDef_HEAD_INIT,
     "gamred_native",
-    "Python interface for the MATLAB fetch_thresholds function",
+    "Python interface for the MATLAB fetch_thresholds & find_gaussian_mixtures functions",
     -1,
     GamredNativeMethods
 };