add em and draft nll

bsi_zoo/estimators.py

@@ -667,3 +667,77 @@ def champagne(L, y, cov=1.0, alpha=0.2, max_iter=1000, max_iter_reweighting=10):
     x[active_set, :] = x_bar
 
     return x
+
+
+def lemur(L, y, max_iter=1000, max_iter_em=10):
+    """Latent EM Unsupervised Regression based on https://ieeexplore.ieee.org/document/9746697
+
+    Parameters
+    ----------
+    L : array, shape (n_sensors, n_sources)
+        lead field matrix modeling the forward operator or dictionary matrix
+    y : array, shape (n_sensors,)
+        measurement vector, capturing sensor measurements
+    max_iter : int, optional
+        The maximum number of inner loop iterations
+    max_iter_reweighting : int, optional
+        Maximum number of reweighting steps i.e outer loop iterations
+
+    Returns
+    -------
+    x : array, shape (n_sources,)
+        Parameter vector, e.g., source vector in the context of BSI (x in the cost
+        function formula).
+
+    References
+    ----------
+    XXX
+    """
+    n_sensors, n_sources = L.shape
+    _, n_times = y.shape
+
+    def moments(Y):
+        """Moments identification method for Gaussian mixture."""
+
+        m2, m4, m6 = np.mean(Y ** 2), np.mean(Y ** 4) / 3, np.mean(Y ** 6) / 15
+
+        A, B = m4 - m2 ** 2, m6 / m2 - m2 ** 2
+        C = (B / A - 3) * m2
+
+        D = (C ** 2) / 4 + A
+        D = max(0, D)
+        X = -C / 2 + np.sqrt(D)
+
+        return (X ** 2 / (A + X ** 2),A / X + X, abs(m2 - X))
+
+    def em_step(obs, param):
+        """EM update with x as complete data."""
+
+        rho = (1-param[0])/param[0]
+        mu = param[1]/(param[1]+param[2])
+
+        s_x = param[1]**2
+        s_b = param[2]**2
+
+        phi_k = 1/(
+            1 + rho*np.sqrt(param[1]/param[2] + 1) * np.exp( -np.sum(obs**2,axis=1)/2 *mu/param[2] )
+            )
+
+        p = np.mean(phi_k)
+        s_x = mu*param[2] + mu**2/p * np.mean(phi_k*np.mean(obs**2,axis=1) )
+        s_b = np.mean(obs**2) - 2 * mu * np.mean(phi_k*np.mean(obs**2,axis=1)) + p*s_x**2
+
+        X_eap = (phi_k * obs.T).T * s_x / (s_x + s_b)
+
+        return ([p, s_x, s_b], X_eap,phi_k)
+
+    x = L.T@y # initialisation of X
+    theta_p = [0,0,0] # initialisation of theta        
+    norm = np.linalg.norm([email protected],2)
+
+    for k_inner in range(max_iter):
+        z = x + L.T@(y - L@x)/norm# Gradient descent
+        theta = moments(z)# Initialisation of the EM (feel free to find better ones !)
+        for k_em in range(max_iter_em):
+            theta, u, phi = em_step(z, theta)# EM updates
+    return x

bsi_zoo/metrics.py

@@ -80,3 +80,23 @@ def euclidean_distance(x, x_hat, *args, **kwargs):
     )
 
     return np.mean(euclidean_distance)
+
+
+def nll(x, x_hat, *args, **kwargs):
+    y = kwargs["y"]
+    L = kwargs["L"]
+    cov = kwargs["cov"]
+    orientation_type = kwargs["orientation_type"]
+    subject = kwargs["subject"]
+    nnz = kwargs["nnz"]
+
+    active_set = _get_active_nnz(x, x_hat, orientation_type, subject, nnz)#kwargs["active_set"]
+
+    # Marginal NegLogLikelihood score upon estimation of the support:
+    # ||(cov + L Q L.T)^-1/2 y||^2_F  + log|cov + L Q L.T| with Q the support matrix
+    q = np.zeros(x.shape[0])
+    q[active_set] = 1
+    Q = np.diag(q)
+    cov_y = cov + L@[email protected]
+    # To take into account the knowledge on nnz you need to add +2log((n_sources-nnz)/nnz)||q||_0
+    return np.linalg.norm(np.linalg.sqrt(np.linalg.inv(cov_y)),ord='fro')**2 + np.log(np.linalg.det(cov_y))