Obliczanie nieznanej wartości p

Niedawno debugowałem skrypt R i znalazłem coś bardzo dziwnego, autor zdefiniował własną funkcję wartości p

pval <- function(x, y){
    if (x+y<20) { # x + y is small, requires R.basic
        p1<- nChooseK(x+y,x) * 2^-(x+y+1);
        p2<- nChooseK(x+y,y) * 2^-(x+y+1);
        pvalue = max(p1, p2)
    }
    else { # if x+y is large, use approximation
        log_p1 <- (x+y)*log(x+y) - x*log(x) - y*log(y) - (x+y+1)*log(2);
        pvalue<-exp(log_p1);
    }
    return(pvalue)
}

Gdzie X i Y są dodatnimi wartościami większymi niż 0. Przypadek <20 wydaje się być obliczeniem pewnego rodzaju rozkładu hipergeometrycznego (coś podobnego do testu Fishera?) I czy ktoś wie, jakie jest inne obliczenie? Jako sidenote, próbuję zoptymalizować ten kod, więc próbuję znaleźć odpowiednią funkcję R do wywołania i zastąpienia go.

Edycja: Wzór wyszczególniania papieru do obliczania wartości p można znaleźć tutaj (należy kliknąć pdf, aby zobaczyć wzory) Metody zaczynają się na stronie 8 pdf, a wzór można znaleźć na stronie 9 w punkcie (1). Zakładają, że rozkład jest Poissonem.

Druga rzecz wygląda na to, że jest to przybliżenie obliczeń zastosowanych w x+y < 20przypadku, ale oparte na przybliżeniu Stirlinga .

Zwykle, gdy jest on używany do tego rodzaju przybliżenia, ludzie używają co najmniej następnego dodatkowego terminu (współczynnik w przybliżeniu dla ), Co poprawiłoby względne przybliżenie znacznie dla małego .$\sqrt{2\pi n}$$n!$$n$

Na przykład, jeżeli i są równe 10, pierwsze obliczenie daje około 0,088 zaś zbliżanie gdy współczynnik znajduje się we wszystkich warunkach wynosi około 0,089, na tyle blisko, dla większości celów ... ale pominięcie tego terminu w przybliżeniu daje 0,5 - co tak naprawdę nie jest wystarczająco blisko! Autor tej funkcji najwyraźniej nie zadał sobie trudu, aby sprawdzić dokładność swojego przybliżenia w przypadku granicy.$x$$y$$\sqrt{2\pi n}$

W tym celu autor powinien prawdopodobnie po prostu wywołać funkcję wbudowaną lgamma- w szczególności, używając tego zamiast tego, co ma do log_p1:

log_p1 <- lgamma(x+y+1)-lgamma(x+1)-lgamma(y+1)-(x+y+1)*log(2)

co skutkuje odpowiedzią, którą próbuje przybliżyć (ponieważ lgamma(x+1)faktycznie zwraca , tą właśnie rzeczą stara się przybliżyć - słabo - poprzez przybliżenie Stirlinga).$\log(x!)$

Podobnie, nie jestem pewien, dlaczego autor nie korzysta z wbudowanej choosefunkcji w pierwszej części, która wchodzi w standardowy rozkład R. Jeśli o to chodzi, odpowiednia funkcja dystrybucji jest prawdopodobnie również wbudowana.

Tak naprawdę nie potrzebujesz dwóch osobnych przypadków; ten lgammadziała dobrze aż do najmniejszych wartości. Z drugiej strony choosefunkcja działa dla całkiem dużych wartości (np. choose(1000,500)Działa dobrze). Bezpieczniejszym rozwiązaniem jest prawdopodobnie lgamma, chociaż trzeba by mieć dość dużą i , zanim było to problemem.$x$$y$

Przy większej ilości informacji powinna istnieć możliwość zidentyfikowania źródła testu. Domyślam się, że autor skądś go wziął, więc powinno być możliwe wyśledzenie go. Czy masz na to jakiś kontekst?

Kiedy mówisz „optymalizuj”, masz na myśli uczynienie go szybszym, krótszym, łatwiejszym w utrzymaniu lub czymś innym?

Edytuj po szybkim przeczytaniu na papierze:

Autorzy wydają się mylnie pod wieloma względami. Dokładny test Fishera nie zakładamy, że marginesy są stałe, to po prostu warunki na nich, który nie jest to samo, w ogóle, jak to opisano, na przykład, tutaj , z odniesieniami. Rzeczywiście, wydają się oni prawie zupełnie nieświadomi debaty na temat uzależnienia od marż i dlaczego tak się dzieje. Linki tam są warte przeczytania.

[Przechodzą od „testu Fishera zawsze bardziej konserwatywnego niż nasz” do twierdzenia, że ​​test Fishera jest zbyt konserwatywny ... co niekoniecznie następuje, chyba że złym jest warunek . Musieliby to ustalić, ale biorąc pod uwagę, że statystycy o to kłócą się od około 80 lat, a ci autorzy wydają się nie wiedzieć, dlaczego uwarunkowanie jest zrobione, nie sądzę, że ci faceci doszli do sedna tego problemu .]

Autorzy artykułu zdają się przynajmniej rozumieć, że podane przez nich prawdopodobieństwa muszą się kumulować, aby uzyskać wartości p; na przykład w pobliżu środka pierwszej kolumny strony 5 (moje podkreślenie):

Istotność statystyczna według dokładnego testu Fishera dla takiego wyniku wynosi 4,6% (wartość P dla dwóch ogonów, tj. Prawdopodobieństwo wystąpienia takiej tabeli w hipotezie, że częstotliwości aktyny EST są niezależne od bibliotek cDNA). Dla porównania wartość P obliczona z postaci skumulowanej (równanie 9, patrz Metody) równania 2 (tj. Dla względnej częstotliwości aktynowych EST w obu bibliotekach, biorąc pod uwagę, że co najmniej 11 pokrewnych EST jest obserwowanych w biblioteka wątroby po dwóch zaobserwowano w bibliotece mózgu) wynosi 1,6%.

(choć nie jestem pewien, czy zgadzam się z ich obliczeniem wartości tam; musiałbym dokładnie sprawdzić, co faktycznie robią z drugim ogonem).

Nie sądzę, że program to robi.

Uważaj jednak, że ich analiza nie jest standardowym testem dwumianowym; używają argumentu Bayesa, aby uzyskać wartość p w teście, który często bywa. Wydaje mi się również - nieco dziwnie, moim zdaniem - warunkować , a nie . Oznacza to, że muszą kończyć się czymś w rodzaju ujemnego dwumianu zamiast dwumianu, ale uważam, że gazeta jest naprawdę źle zorganizowana i strasznie źle wyjaśniona (i jestem przyzwyczajony do opracowywania tego, co dzieje się w dokumentach statystycznych), więc jestem nie będę pewien, dopóki nie przejdę ostrożnie.$x$$x+y$

Nie jestem nawet przekonany, że suma ich prawdopodobieństw wynosi w tym momencie 1.

Jest tu o wiele więcej do powiedzenia, ale pytanie nie dotyczy papieru, ale implementacji w programie.

-

W każdym razie wynik jest taki, że przynajmniej papier poprawnie identyfikuje, że wartości p składają się z sumy prawdopodobieństw takich jak te w równaniu 2, ale program tego nie robi . (Patrz eqn 9a i 9b w sekcji Metody w dokumencie).

Kod jest po prostu zły.

[Mógłbyś użyć pbinom, jak sugerowałby komentarz @ Whubera, do obliczenia poszczególnych prawdopodobieństw (ale nie ogona, ponieważ nie jest to test dwumianowy, ponieważ je konstruują), ale wtedy jest dodatkowy czynnik 1/2 w ich równaniu 2, więc jeśli chcesz powielić wyniki w pracy, musisz je zmienić.]

Możesz go uzyskać, z pewnymi skrzypkami, od pnbinom-

Zwyczajowe negatywu dwumianowego są albo liczbę prób do sukces lub liczbę niepowodzeń do sukcesem. Te dwa są równoważne; Wikipedia podaje tutaj drugą formę . Funkcja prawdopodobieństwa to:$k^\rm{th}$$k^\rm{th}$

Równanie 2 na p4 (a więc także równanie 1 na p3) jest dwumianem ujemnym, ale przesuniętym o 1. Niech , i .$p = N_1/(N_1+N_2)$$k=x$$r = y+1$

Niepokoi mnie to, że ponieważ limity nie zostały podobnie przesunięte, ich prawdopodobieństwa mogą nawet nie wzrosnąć do 1.$y$

To by było złe.