diff --git a/corex.py b/corex.py
index 824be8c..606ebad 100644
--- a/corex.py
+++ b/corex.py
@@ -494,7 +494,7 @@ def marginal_p(self, xi, thetai):
         if self.marginal_description == 'gaussian':
             mu, sig = thetai  # mu, sig have size m by k
             xi = xi.reshape((-1, 1, 1))
-            return (-(xi - mu)**2 / (2. * sig) - 0.5 * np.log(2 * np.pi * sig)).transpose((1, 0, 2))  # log p(xi|yj)
+            return (-(xi.data - mu)**2 / (2. * sig) - 0.5 * np.log(2 * np.pi * sig)).transpose((1, 0, 2))  # log p(xi|yj) # intentionally removing mask on xi prior to calculation to avoid bug: https://github.com/numpy/numpy/issues/18744
 
         elif self.marginal_description == 'discrete':
             # Discrete data: should be non-negative integers starting at 0: 0,...k. k < 32 because of np.choose limits
@@ -510,7 +510,7 @@ def estimate_parameters(self, xi, p_y_given_x):
             n_obs = np.sum(p_y_given_x, axis=1).clip(0.1)  # m, k
             mean_ml = np.einsum('i,jik->jk', xi, p_y_given_x, optimize=False) / n_obs  # ML estimate of mean of Xi
             sig_ml = np.einsum('jik,jik->jk', (xi[np.newaxis, :, np.newaxis] - mean_ml[:, np.newaxis, :])**2, p_y_given_x, optimize=False) / (n_obs - 1).clip(0.01)  # UB estimate of sigma^2(variance)
-            sig_ml = sig_ml.clip(0.25) # To prevent divide by zero in marginal_p function
+            sig_ml = sig_ml.clip(1e-200) # To prevent divide by zero in marginal_p function
             
             if not self.smooth_marginals:
                 return np.array([mean_ml, sig_ml])  # FOR EACH Y_j = k !!