Cómo explicar la división de una distribución bimodal basada en la estimación de la densidad del núcleo

Question

Tengo un conjunto de datos de población bimodal. Contiene un pico más pequeño, que se considera "malo", y un pico más grande. Intento separar la parte mala de los datos del resto de los datos. Lo que hice fue: primero hice una estimación de la densidad del núcleo, luego encontré el máximo local de este pequeño pico, y el mínimo local de la fosa entre dos picos, luego tomé el punto medio (media aritmética de las coordenadas x) de ellos, y lo definí como un corte. Todo lo que está por debajo de este límite se considera "malo". La razón por la que tomé el punto medio en lugar de la fosa es porque intenté ser más conservador.

Ahora me gustaría preguntar: ¿Es razonable lo que hice? Si la respuesta es afirmativa, ¿cómo puedo explicar mi acción de una manera que favorezca a los estadísticos? Si no, ¿cómo puedo cambiar? (Cualquier otro método es bienvenido, especialmente los implementados en R.) ¡Gracias!

Esta es la cifra.

Answer 1

Se podría ajustar un modelo de mezcla de dos componentes utilizando http://cran.r-project.org/web/packages/mixtools/index.html . Prueba a utilizar normalmixEM. A continuación, podría seguir las sugerencias de Erich Schubert y encontrar la región en la que Pr[punto de datos se generó a partir del componente con la media más pequeña] >= 0,50.

Edición: ejemplo de código R:

library(mixtools)

simulate <- function(lambda=0.3, mu=c(0, 4), sd=c(1, 1), n.obs=10^5) {
    x1 <- rnorm(n.obs, mu[1], sd[1])
    x2 <- rnorm(n.obs, mu[2], sd[2])    
    return(ifelse(runif(n.obs) < lambda, x1, x2))
}

x <- simulate()

model <- normalmixEM(x=x, k=2)
index.lower <- which.min(model$mu)  # Index of component with lower mean

find.cutoff <- function(proba=0.5, i=index.lower) {
    ## Cutoff such that Pr[drawn from bad component] == proba
    f <- function(x) {
        proba - (model$lambda[i]*dnorm(x, model$mu[i], model$sigma[i]) /
                     (model$lambda[1]*dnorm(x, model$mu[1], model$sigma[1]) + model$lambda[2]*dnorm(x, model$mu[2], model$sigma[2])))
        }
        return(uniroot(f=f, lower=-10, upper=10)$root)  # Careful with division by zero if changing lower and upper
}

cutoffs <- c(find.cutoff(proba=0.5), find.cutoff(proba=0.75))  # Around c(1.8, 1.5)

hist(x)
abline(v=cutoffs, col=c("red", "blue"), lty=2)

Answer 2

Probablemente tendría más sentido si también se estimara la "altura" (en realidad, más bien el peso) de ambos, y luego se fijara el umbral en el punto de inflexión.

Es decir, modelar los datos como $$p_1 \cdot pdf(x, \mu_1, \sigma_1) + p_2 \cdot pdf(x, \mu_2, \sigma_2)$$

y establecer el umbral en $x$ donde $$p_1 \cdot pdf(x, \mu_1, \sigma_1) = p_2 \cdot pdf(x, \mu_2, \sigma_2)$$

es decir, el objeto tiene la misma probabilidad de pertenecer a ambas clases.

Todavía puedes añadir un parámetro para afinar lo conservador que es tu método, por ejemplo, usando $$p_1 \cdot pdf(x, \mu_1, \sigma_1) = c\cdot p_2 \cdot pdf(x, \mu_2, \sigma_2)$$

donde $c=2$ pondría el doble de peso en la segunda distribución.

Answer 3

Estoy utilizando este ejemplo y a veces me sale este error

Error en uniroot(f = f, inferior = -10, superior = 10) : Los valores de f() en los puntos finales no son de signo contrario

Así que cambié el valor inferior a -1 y para algunos de los conjuntos de datos se arregló, pero sigue dando errores en otros. No estoy seguro de si se puede establecer dinámicamente basado en el vector de entrada (es decir, x)?