Accounting for noise in 2D Gaussian model

Question

I need to fit a 2D gaussian embedded into substantial uniform noise, as shown in the left plot below. I tried using sklearn.mixture.GaussianMixture with two components (code at the bottom), but this obviously fails as shown in the right plot below.

I want to assign probabilities to each element of belonging to the 2D Gaussian and to the uniform background noise. This seems like a simple enough task but I've found no "simple" way to do it.

Any advices? It doesn't need to be GMM, I'm open to other methods/packages.

---

import numpy as np
import matplotlib.pyplot as plt
from sklearn import mixture

# Generate 2D Gaussian data
N_c = 100
xy_c = np.random.normal((.5, .5), .05, (N_c, 2))

# Generate uniform noise
N_n = 1000
xy_n = np.random.uniform(.0, 1., (N_n, 2))

# Combine into a single data set
data = np.concatenate([xy_c, xy_n])

# fit a Gaussian Mixture Model with two components
model = mixture.GaussianMixture(n_components=2, covariance_type='full')
model.fit(data)
probs = model.predict_proba(data)
labels = model.predict(data)
# Separate the two clusters for plotting
msk0 = labels == 0
c0, p0 = data[msk0], probs[msk0].T[0]
msk1 = labels == 1
c1, p1 = data[msk1], probs[msk1].T[1]

# Plot
plt.subplot(121)
plt.scatter(*xy_n.T, c='b', alpha=.5)
plt.scatter(*xy_c.T, c='r', alpha=.5)
plt.xlim(0., 1.)
plt.ylim(0., 1.)

plt.subplot(122)
plt.scatter(*c0.T, c=p0, alpha=.75)
plt.scatter(*c1.T, c=p1, alpha=.75)
plt.colorbar()
# display predicted scores by the model as a contour plot
X, Y = np.meshgrid(np.linspace(0., 1.), np.linspace(0., 1.))
XX = np.array([X.ravel(), Y.ravel()]).T
Z = -model.score_samples(XX)
Z = Z.reshape(X.shape)
plt.contour(X, Y, Z)

plt.show()

Answer 1

I think kernel density can help you to localize the gaussian and exclude point outside of it (e.g in area with lesser densities)

Here is an example code :

import numpy as np
import matplotlib.pyplot as plt
from sklearn import mixture
from sklearn.neighbors import KernelDensity


# Generate 2D Gaussian data
N_c = 100
xy_c = np.random.normal((.2, .2), .05, (N_c, 2))

# Generate uniform noise
N_n = 1000
xy_n = np.random.uniform(.0, 1., (N_n, 2))

# Combine into a single data set
data = np.concatenate([xy_c, xy_n])
print(data.shape)

model = KernelDensity(kernel='gaussian',bandwidth=0.05)
model.fit(data)
probs = model.score_samples(data)

# Plot
plt.subplot(131)
plt.scatter(*xy_n.T, c='b', alpha=.5)
plt.scatter(*xy_c.T, c='r', alpha=.5)

plt.xlim(0., 1.)
plt.ylim(0., 1.)

# plot kernel score
plt.subplot(132)
plt.scatter(*data.T, c=probs, alpha=.5)

# display predicted scores by the model as a contour plot
X, Y = np.meshgrid(np.linspace(0., 1.), np.linspace(0., 1.))
XX = np.array([X.ravel(), Y.ravel()]).T
Z = -model.score_samples(XX)
Z = Z.reshape(X.shape)
plt.contour(X, Y, Z)
plt.xlim(0,1)
plt.ylim(0,1)

# plot kernel score with threshold
plt.subplot(133)
plt.scatter(*data.T, c=probs>0.5, alpha=.5) # here you can adjust the threshold
plt.colorbar()
plt.xlim(0,1)
plt.ylim(0,1)

And this is the output figure :

I changed the center of the gaussian to ensure my code was working. The right panel display the kernel score with a threshold, which can be use in your case to filter out the noisy data outside of the gaussian, but you can't filter the noise inside the gaussian.