Dobroć dopasowania do histogramów 2D

Mam dwa zestawy danych reprezentujących parametry gwiazd: obserwowany i modelowany. Za pomocą tych zestawów tworzę tak zwany schemat dwukolorowy (TCD). Próbkę można zobaczyć tutaj:

Być obserwowane dane i B dane wydobyte z modelu (nieważne czarne linie, kropki reprezentują dane) Mam tylko jedno A schemat, ale może produkować tyle różnych B schematów, jak chcę, i to, co jest mi potrzebne aby utrzymać ten, który najlepiej pasuje A .

Potrzebuję więc niezawodnego sposobu sprawdzenia poprawności dopasowania schematu B (model) do schematu A (zaobserwowany).

Teraz tworzę histogram 2D ​​lub siatkę (tak to nazywam, może ma bardziej odpowiednią nazwę) dla każdego diagramu, dzieląc obie osie (100 pojemników na każdą) Następnie przechodzę przez każdą komórkę siatki i znajduję absolutną różnicę w liczbie między A i B dla tej konkretnej komórki. Po przejściu przez wszystkie komórki, to suma wartości dla każdej komórki, a więc kończy się z jednej pozytywnej parametru reprezentującego dobroć dopasowania ( ) pomiędzy A i B . Im najbliżej zera, tym lepsze dopasowanie. Zasadniczo tak wygląda ten parametr:$gf$

a i j i j b i $gf = \sum_{ij} |a_{ij}-b_{ij}|$; gdzie jest liczba gwiazd na schemacie A dla danej komórki (oznaczonej przez ) oraz jest numerem B .$a_{ij}$$ij$$b_{ij}$

Tak wyglądają te różnice w liczbie w każdej komórce w utworzonej przeze mnie siatce (zwróć uwagę, że nie używam wartości bezwzględnych w tym obrazie, ale zrobić z nich korzystać przy obliczaniu parametr):( a i j - b i j ) g $(a_{ij}-b{ij})$$(a_{ij}-b{ij})$$gf$

Problem polega na tym, że doradzono mi, że może to nie być dobry estymator, głównie dlatego, że poza tym, że dopasowanie jest lepsze niż inne, ponieważ parametr jest niższy , naprawdę nie mogę nic więcej powiedzieć.

Ważne :

(dzięki @PeterEllis za podniesienie tego)

1- Punkty w B nie są związane jeden do jednego z punktów A . To ważna rzecz, aby pamiętać przy poszukiwaniu najlepszego dopasowania: liczba punktów A i B to nie koniecznie taka sama, a dobroć dopasowania testu należy również uwzględnić tej rozbieżności i spróbować go zminimalizować.

2- Liczba punktów w każdym zestawie danych B (model wyjściowy), który próbuję dopasować do A, nie jest stała.

W niektórych przypadkach widziałem test Chi-Squared :

O i E $\sum_i (O_i-E_i)^2/E_i$ ; gdzie to częstotliwość obserwowana (model), a to częstotliwość oczekiwana (obserwacja).$O_i$$E_i$

ale problem polega na tym: co mam zrobić, jeśli wynosi zero? Jak widać na powyższym obrazku, jeśli utworzę siatkę tych diagramów w tym zakresie, będzie wiele komórek, w których wynosi zero.E $E_i$$E_i$

Przeczytałem również, że niektórzy ludzie zalecają test prawdopodobieństwa logarytmicznego Poissona do zastosowania w przypadkach takich jak ten, w których uczestniczą histogramy. Jeśli to prawda, naprawdę doceniłbym to, gdyby ktoś mógł poinstruować mnie, jak korzystać z tego testu w tym konkretnym przypadku (pamiętaj, że moja wiedza na temat statystyk jest dość słaba, więc proszę, zachowaj to tak proste, jak to możliwe :)

OK, gruntownie poprawiłem tę odpowiedź. Wydaje mi się, że zamiast binować dane i porównywać liczby w każdym bin, sugestia, którą zakopałem w mojej oryginalnej odpowiedzi na temat dopasowania oszacowania gęstości jądra 2d i porównania ich, jest znacznie lepszym pomysłem. Co więcej, w pakiecie ks Tarna Duonga dla R znajduje się funkcja kde.test (), która robi to łatwo.

Sprawdź dokumentację kde.test, aby uzyskać więcej informacji i argumentów, które możesz poprawić. Ale w zasadzie robi dokładnie to, czego chcesz. Zwracana wartość p to prawdopodobieństwo wygenerowania dwóch zestawów danych, które porównujesz w ramach hipotezy zerowej, że były one generowane z tego samego rozkładu. Więc im wyższa wartość p, tym lepsze dopasowanie między A i B. Zobacz mój przykład poniżej, gdzie łatwo zauważyć, że B1 i A są różne, ale że B2 i A są prawdopodobnie takie same (czyli w jaki sposób zostały wygenerowane) .

# 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

MOJA ORYGINALNA ODPOWIEDŹ PONIŻEJ JEST UTRZYMYWANA TYLKO PONIEWAŻ, ŻE SĄ TERAZ ŁĄCZA Z NIMI POZOSTAŁE, KTÓRE NIE UCZYNIĄ

Po pierwsze, mogą istnieć inne sposoby rozwiązania tego problemu.

Justel i wsp. Przedstawili wielowymiarowe rozszerzenie testu dobroci dopasowania Kołmogorowa-Smirnowa, które moim zdaniem można by zastosować w twoim przypadku, aby sprawdzić, jak dobrze każdy zestaw modelowanych danych pasuje do oryginału. Nie mogłem znaleźć implementacji tego (np. W R), ale może nie wyglądałem wystarczająco mocno.

Alternatywnie może istnieć sposób, aby to zrobić, dopasowując kopułę zarówno do oryginalnych danych, jak i do każdego zestawu modelowanych danych, a następnie porównując te modele. Istnieją implementacje tego podejścia w R i innych miejscach, ale nie jestem specjalnie z nimi zaznajomiony, więc nie próbowałem.

Ale aby odpowiedzieć bezpośrednio na twoje pytanie, przyjęte podejście jest rozsądne. Sugeruje się kilka punktów:

• Chyba że twój zestaw danych jest większy niż się wydaje, myślę, że siatka 100 x 100 to zbyt wiele pojemników. Intuicyjnie mogę sobie wyobrazić, że wyciąganie wniosków z różnych zestawów danych jest bardziej odmienne niż tylko z powodu precyzji twoich pojemników, co oznacza, że ​​masz wiele pojemników z małą liczbą punktów, nawet gdy gęstość danych jest wysoka. Jednak ostatecznie jest to kwestia osądu. Z pewnością sprawdziłbym twoje wyniki za pomocą różnych podejść do binowania.

• Po zakończeniu binowania i przekonwertowaniu danych na (w efekcie) tabelę zdarzeń z dwiema kolumnami i liczbą wierszy równą liczbie pojemników (10 000 w twoim przypadku), masz standardowy problem z porównaniem dwóch kolumn liczy. Sprawdziłby się albo test chi-kwadrat, albo dopasowanie jakiegoś modelu Poissona, ale jak mówisz, istnieje niezręczność z powodu dużej liczby zer. Każdy z tych modeli jest zwykle dopasowany, minimalizując sumę kwadratów różnicy, ważoną odwrotnością oczekiwanej liczby zliczeń; gdy zbliża się do zera, może powodować problemy.

Edytuj - reszty tej odpowiedzi nie uważam już za właściwe podejście.

Myślę, że dokładny test Fishera może nie być przydatny ani odpowiedni w tej sytuacji, w której krańcowe sumy wierszy w tabeli krzyżowej nie są ustalone. Daje wiarygodną odpowiedź, ale trudno mi pogodzić jej użycie z oryginalną pochodną z projektu eksperymentalnego. Zostawiam tutaj oryginalną odpowiedź, więc komentarze i pytania uzupełniające mają sens. Ponadto nadal może istnieć sposób odpowiedzi na pożądane podejście OP polegające na binowaniu danych i porównywaniu pojemników za pomocą jakiegoś testu opartego na średnich różnicach bezwzględnych lub kwadratowych. W takim podejściu nadal wykorzystano by , o której mowa poniżej, i sprawdzono niezależność, tj. Szukając wyniku, w którym kolumna A miałaby takie same proporcje jak kolumna B.$n_g\times 2$

Podejrzewam, że rozwiązaniem powyższego problemu byłoby użycie dokładnego testu Fishera , stosując go do , gdzie jest całkowitą liczbą pojemników. Chociaż pełne obliczenia prawdopodobnie nie będą praktyczne ze względu na liczbę wierszy w tabeli, można uzyskać dobre oszacowania wartości p za pomocą symulacji Monte Carlo (implementacja R testu Fishera daje tę opcję jako opcję dla tabel, które są większy niż 2 x 2 i podejrzewam, że robią to inne pakiety). Te wartości p są prawdopodobieństwem, że drugi zestaw danych (z jednego z twoich modeli) ma taki sam rozkład w pojemnikach jak oryginał. Dlatego im wyższa wartość p, tym lepsze dopasowanie. $n_g \times 2$$n_g$

Symulowałem niektóre dane, aby wyglądać trochę jak twoje, i odkryłem, że to podejście było dość skuteczne w identyfikacji, które z moich zestawów danych „B” zostały wygenerowane z tego samego procesu co „A”, a które nieco się różniły. Z pewnością bardziej skuteczny niż gołym okiem.

• Przy takim podejściu do testowania niezależności zmiennych w tabeli kontyngencji nie ma znaczenia, że ​​liczba punktów w A jest różna od liczby w B (chociaż zauważ, że jest to$n_g \times 2$problem, jeśli użyjesz tylko sumy różnic bezwzględnych lub różnic kwadratowych, jak pierwotnie proponujesz). Jednak nie ma znaczenia, że ​​każda z twoich wersji B ma inną liczbę punktów. Zasadniczo większe zestawy danych B będą miały tendencję do zwracania niższych wartości p. Mogę wymyślić kilka możliwych rozwiązań tego problemu. 1. Mógłbyś zredukować wszystkie swoje zestawy B danych do tego samego rozmiaru (rozmiar najmniejszego ze swoich zestawów B), pobierając losową próbkę tego rozmiaru ze wszystkich zestawów B, które są większe niż ten rozmiar. 2. Możesz najpierw dopasować dwuwymiarowe oszacowanie gęstości jądra do każdego ze swoich zestawów B, a następnie zasymulować dane z tego oszacowania o równych rozmiarach. 3. możesz użyć jakiejś symulacji, aby obliczyć stosunek wartości p do wielkości i użyć tego do „poprawienia” wartości p otrzymane z powyższej procedury, więc są one porównywalne. Prawdopodobnie istnieją też inne alternatywy. To, co zrobisz, będzie zależeć od tego, jak wygenerowano dane B, jak różne są rozmiary itp.

Mam nadzieję, że to pomaga.