Prueba de bondad de ajuste para una mezcla en R

Question

Acabo de estimar los parámetros de una mezcla de dos gaussianos con diferentes medias y diferentes sigmas, me gustaría probar si los datos se ajustan bien a la forma explícita de la mezcla, ¿necesito necesariamente simular los datos de la mezcla o cómo puedo probar la bondad del ajuste? Estoy utilizando el paquete mixtools.

Answer 1

Se puede escribir una función que calcule el valor relevante para una prueba dada bajo la hipótesis nula basándose en las salidas de, por ejemplo, normalmixEM y luego hacer la prueba con eso. Por ejemplo, para la prueba de Kolmogorov-Smirnov, necesitamos la FCD dado un conjunto de parámetros:

# CDF of mixture of two normals
pmnorm <- function(x, mu, sigma, pmix) {
  pmix[1]*pnorm(x,mu[1],sigma[1]) + (1-pmix[1])*pnorm(x,mu[2],sigma[2])
}

A continuación, realizamos la prueba K-S de la forma habitual:

# Sample run
x <- c(rnorm(50), rnorm(50,2))

foo <- normalmixEM(x)
test <- ks.test(x, pmnorm, mu=foo$mu, sigma=foo$sigma, pmix=foo$lambda)
test

    One-sample Kolmogorov-Smirnov test

data:  x 
D = 0.0559, p-value = 0.914
alternative hypothesis: two-sided 

Siempre teniendo en cuenta el hecho de que estimamos los parámetros a partir de los mismos datos que estamos utilizando para hacer la prueba, sesgando así la prueba hacia el no rechazo de H0.

Podemos superar este último sesgo hasta cierto punto mediante un bootstrap paramétrico, generando muchas muestras a partir de una mezcla de normales parametrizadas por las estimaciones de normalmixEM y, a continuación, estimar los parámetros de las muestras y calcular los estadísticos de la prueba para cada muestra utilizando los parámetros estimados. Bajo esta construcción, la hipótesis nula es siempre verdadera. En el código de abajo, estoy ayudando al algoritmo EM empezando por los parámetros verdaderos de la muestra, lo que también es hacer un poco de trampa, ya que hace que el algoritmo EM tenga más probabilidades de encontrar valores cercanos a los verdaderos que con la muestra original, pero reduce enormemente el número de mensajes de error.

# Bootstrap estimation of ks statistic distribution
N <- length(x)
ks.boot <- rep(0,1000)
for (i in 1:1000) {
  z <- rbinom(N, 1, foo$lambda[1])
  x.b <- z*rnorm(N, foo$mu[1], foo$sigma[1]) + (1-z)*rnorm(N, foo$mu[2], foo$sigma[2])
  foo.b <- normalmixEM(x.b, maxit=10000, lambda=foo$lambda, mu=foo$mu, sigma=foo$sigma)
  ks.boot[i] <- ks.test(x.b, pmnorm, mu=foo.b$mu, sigma=foo.b$sigma, pmix=foo.b$lambda)$statistic
}

mean(test$statistic <= ks.boot)
[1] 0.323

Así, en lugar de un valor p de 0,914, obtenemos un valor p de 0,323. Interesante, pero no especialmente importante en este caso.