Bondad de ajuste para histogramas 2D

Tengo dos conjuntos de datos que representan los parámetros de las estrellas: uno observado y uno modelado. Con estos conjuntos creo lo que se llama un diagrama de dos colores (TCD). Una muestra se puede ver aquí:

A son los datos observados y B los datos extraídos del modelo (no importa las líneas negras, los puntos representan los datos) Solo tengo un diagrama A , pero puedo producir tantos diagramas B como quiera, y lo que necesito es para mantener el que mejor se adapte a .

Entonces, lo que necesito es una forma confiable de verificar la bondad de ajuste del diagrama B (modelo) al diagrama A (observado).

En este momento, lo que hago es crear un histograma o cuadrícula 2D (así es como lo llamo, tal vez tiene un nombre más propio) para cada diagrama agrupando ambos ejes (100 contenedores para cada uno) Luego paso por cada celda de la cuadrícula y encuentro la diferencia absoluta en los recuentos entre A y B para esa celda en particular. Después de haber pasado por todas las celdas, sumo los valores para cada celda y así termino con un solo parámetro positivo que representa la bondad de ajuste ( ) entre A y$gf$ B . Cuanto más cercano a cero, mejor será el ajuste. Básicamente, así es como se ve ese parámetro:

; donde un i j es el número de estrellas en el diagramaApara esa célula en particular (determinado por i j ) y b i j es el número deB.$gf = \sum_{ij} |a_{ij}-b_{ij}|$$a_{ij}$$ij$$b_{ij}$

Así es como se ven esas diferencias en los recuentos en cada celda en la cuadrícula que creo (tenga en cuenta que no estoy usando valores absolutos de ( a i j - b i j ) en esta imagen, pero I hago utilizarlos en el cálculo de la g f parámetro):$(a_{ij}-b{ij})$$(a_{ij}-b{ij})$$gf$

El problema es que me han informado que este podría no ser un buen estimador, principalmente porque aparte de decir que este ajuste es mejor que este otro porque el parámetro es más bajo , realmente no puedo decir nada más.

Importante :

(gracias @PeterEllis por mencionar esto)

1- Puntos en B no están relacionados uno-a-uno con puntos en A . Eso es algo importante a tener en cuenta al buscar el mejor ajuste: el número de puntos en A y B no es necesariamente el mismo y la prueba de bondad de ajuste también debe tener en cuenta esta discrepancia y tratar de minimizarla.

2- El número de puntos en cada conjunto de datos B (salida del modelo) que intento ajustar a A no es fijo.

He visto la prueba Chi-Squared utilizada en algunos casos:

$\sum_i (O_i-E_i)^2/E_i$$O_i$$E_i$

$E_i$$E_i$

Además, he leído que algunas personas recomiendan que se aplique una prueba de Poisson de probabilidad de registro en casos como este donde están involucrados histogramas. Si esto es correcto, realmente agradecería que alguien pudiera instruirme sobre cómo usar esa prueba para este caso en particular (recuerde, mi conocimiento de las estadísticas es bastante pobre, así que manténgalo lo más simple posible :)

OK, he revisado ampliamente esta respuesta. Creo que en lugar de agrupar sus datos y comparar recuentos en cada contenedor, la sugerencia que había enterrado en mi respuesta original de ajustar una estimación de densidad de kernel 2d y compararlos es una idea mucho mejor. Aún mejor, hay una función kde.test () en el paquete ks de Tarn Duong para R que lo hace tan fácil como un pastel.

Consulte la documentación de kde.test para obtener más detalles y los argumentos que puede modificar. Pero básicamente hace exactamente lo que quieres. El valor p que devuelve es la probabilidad de generar los dos conjuntos de datos que está comparando bajo la hipótesis nula de que se estaban generando a partir de la misma distribución. Por lo tanto, cuanto mayor sea el valor p, mejor será el ajuste entre A y B. Vea mi ejemplo a continuación, donde esto fácilmente identifica que B1 y A son diferentes, pero que B2 y A son plausiblemente iguales (que es cómo se generaron) .

# generate some data that at least looks a bit similar
generate <- function(n, displ=1, perturb=1){
    BV <- rnorm(n, 1*displ, 0.4*perturb)
    UB <- -2*displ + BV + exp(rnorm(n,0,.3*perturb))
    data.frame(BV, UB)
}
set.seed(100)
A <- generate(300)
B1 <- generate(500, 0.9, 1.2)
B2 <- generate(100, 1, 1)
AandB <- rbind(A,B1, B2)
AandB$type <- rep(c("A", "B1", "B2"), c(300,500,100))

# plot
p <- ggplot(AandB, aes(x=BV, y=UB)) + facet_grid(~type) + 
    geom_smooth() +     scale_y_reverse() + theme_grey(9)
win.graph(7,3)
p +geom_point(size=.7)

> library(ks)
> kde.test(x1=as.matrix(A), x2=as.matrix(B1))$pvalue
[1] 2.213532e-05
> kde.test(x1=as.matrix(A), x2=as.matrix(B2))$pvalue
[1] 0.5769637

MI RESPUESTA ORIGINAL A CONTINUACIÓN, SE GUARDA SOLO PORQUE AHORA HAY ENLACES DESDE OTRA PARTE QUE NO TENDRÁ SENTIDO

Primero, puede haber otras formas de hacerlo.

Justel y col. han presentado una extensión multivariada de la prueba de bondad de ajuste de Kolmogorov-Smirnov que creo que podría usarse en su caso, para probar qué tan bien se ajusta cada conjunto de datos modelados al original. No pude encontrar una implementación de esto (por ejemplo, en R) pero tal vez no busqué lo suficiente.

Alternativamente, puede haber una manera de hacer esto ajustando una cópula tanto a los datos originales como a cada conjunto de datos modelados, y luego comparando esos modelos. Hay implementaciones de este enfoque en R y otros lugares, pero no estoy especialmente familiarizado con ellos, así que no lo he intentado.

Pero para abordar su pregunta directamente, el enfoque que ha adoptado es razonable. Se sugieren varios puntos:

• A menos que su conjunto de datos sea más grande de lo que parece, creo que una cuadrícula de 100 x 100 es demasiados contenedores. Intuitivamente, puedo imaginar que concluir que los diversos conjuntos de datos son más diferentes de lo que son solo debido a la precisión de sus contenedores significa que tiene muchos contenedores con un bajo número de puntos, incluso cuando la densidad de datos es alta. Sin embargo, esto al final es una cuestión de juicio. Ciertamente verificaría sus resultados con diferentes enfoques para binning.

• Una vez que haya realizado el binning y haya convertido sus datos en (en efecto) una tabla de recuento de contingencia con dos columnas y un número de filas igual al número de bins (10,000 en su caso), tiene un problema estándar de comparar dos columnas de cuentas. Una prueba de Chi cuadrado o la adaptación de algún tipo de modelo de Poisson funcionaría, pero como usted dice, es incómodo debido a la gran cantidad de recuentos cero. Cualquiera de esos modelos normalmente se ajusta minimizando la suma de los cuadrados de la diferencia, ponderada por el inverso del número esperado de conteos; cuando esto se acerca a cero puede causar problemas.

Editar: el resto de esta respuesta ahora ya no creo que sea un enfoque apropiado.

$n_g\times 2$

$n_g \times 2$$n_g$ es el número total de contenedores. Aunque no es probable que un cálculo completo sea práctico debido al número de filas en su tabla, puede obtener buenas estimaciones del valor p utilizando la simulación de Monte Carlo (la implementación R de la prueba de Fisher ofrece esto como una opción para tablas que son más grande que 2 x 2 y sospecho que también lo hacen otros paquetes). Estos valores p son la probabilidad de que el segundo conjunto de datos (de uno de sus modelos) tenga la misma distribución a través de sus bins que el original. Por lo tanto, cuanto mayor sea el valor p, mejor será el ajuste.

Simulé algunos datos para que se parecieran un poco a los suyos y descubrí que este enfoque era bastante efectivo para identificar cuáles de mis conjuntos de datos "B" se generaron a partir del mismo proceso que "A" y cuáles eran ligeramente diferentes. Ciertamente más efectivo que a simple vista.

• $n_g \times 2$ esun problema si usa solo la suma de diferencias absolutas o diferencias cuadradas, como propuso originalmente). Sin embargo, sí importa que cada una de sus versiones de B tenga un número diferente de puntos. Básicamente, los conjuntos de datos B más grandes tenderán a devolver valores p más bajos. Puedo pensar en varias posibles soluciones a este problema. 1. Puede reducir todos sus conjuntos de datos B al mismo tamaño (el tamaño del más pequeño de sus conjuntos B), tomando una muestra aleatoria de ese tamaño de todos los conjuntos B que son más grandes que ese tamaño. 2. Primero puede ajustar una estimación de densidad de núcleo bidimensional a cada uno de sus conjuntos B y luego simular datos de esa estimación que sean de igual tamaño. 3. podría usar algún tipo de simulación para determinar la relación de los valores p con el tamaño y usarla para "corregir" los valores p que obtiene del procedimiento anterior para que sean comparables. Probablemente hay otras alternativas también. El que haga dependerá de cómo se generaron los datos B, cuán diferentes son los tamaños, etc.

Espero que ayude.