Aquí hay un script para usar el modelo de mezcla usando mcluster.
X <- c(rnorm(200, 10, 3), rnorm(200, 25,3), rnorm(200,35,3), rnorm(200,65, 3), rnorm(200,80,5))
Y <- c(rnorm(1000, 30, 2))
plot(X,Y, ylim = c(10, 60), pch = 19, col = "gray40")
require(mclust)
xyMclust <- Mclust(data.frame (X,Y))
plot(xyMclust)
 

En una situación donde hay menos de 5 grupos:
X1 <- c(rnorm(200, 10, 3), rnorm(200, 25,3), rnorm(200,35,3),  rnorm(200,80,5))
Y1 <- c(rnorm(800, 30, 2))
xyMclust <- Mclust(data.frame (X1,Y1))
plot(xyMclust)

 xyMclust4 <- Mclust(data.frame (X1,Y1), G=3)
plot(xyMclust4)

En este caso estamos ajustando 3 grupos. ¿Qué pasa si encajamos 5 grupos?
xyMclust4 <- Mclust(data.frame (X1,Y1), G=5)
plot(xyMclust4)
Se puede forzar a hacer 5 grupos.

También introduzcamos algo de ruido aleatorio: 
X2 <- c(rnorm(200, 10, 3), rnorm(200, 25,3), rnorm(200,35,3),  rnorm(200,80,5), runif(50,1,100 ))
Y2 <- c(rnorm(850, 30, 2))
xyMclust1 <- Mclust(data.frame (X2,Y2))
plot(xyMclust1)
mclustpermite la agrupación basada en modelos con ruido, es decir, observaciones periféricas que no pertenecen a ningún grupo. mclustpermite especificar una distribución previa para regularizar el ajuste a los datos. Se priorControlproporciona una función en mclust para especificar el previo y sus parámetros. Cuando se llama con sus valores predeterminados, invoca otra función llamada defaultPriorque puede servir como plantilla para especificar alternativas anteriores. Para incluir el ruido en el modelado, se debe proporcionar una estimación inicial de las observaciones de ruido a través del componente de ruido del argumento de inicialización en Mclusto mclustBIC.

La otra alternativa sería usar un mixtools paquete que le permita especificar la media y la sigma para cada componente.
X2 <- c(rnorm(200, 10, 3), rnorm(200, 25,3), rnorm(200,35,3),
    rnorm(200,80,5), rpois(50,30))
Y2 <- c(rnorm(800, 30, 2), rpois(50,30))
df <- cbind (X2, Y2)
require(mixtools)
out <- mvnormalmixEM(df, lambda = NULL, mu = NULL, sigma = NULL,
   k = 5,arbmean = TRUE, arbvar = TRUE, epsilon = 1e-08,  maxit = 10000, verb = FALSE)
plot(out, density = TRUE, alpha = c(0.01, 0.05, 0.10, 0.12, 0.15),  marginal = TRUE)

                 
                
Un enfoque estándar son los modelos de mezcla gaussiana que se entrenan mediante el algoritmo EM. Pero como también observa que la cantidad de clústeres puede variar, también puede considerar un modelo no paramétrico como el Dirichlet GMM que también se implementa en scikit-learn.
En R, estos dos paquetes parecen ofrecer lo que necesita,
fuente