Szacowanie rozkładu na podstawie trzech percentyli

23

Jakich metod mogę użyć do wnioskowania o rozkładzie, jeśli znam tylko trzy percentyle?

Na przykład wiem, że w pewnym zbiorze danych piąty percentyl wynosi 8,135, 50 percentyl to 11 259, a 95 percentyl to 23 611. Chcę móc przejść z dowolnej innej liczby do jej percentyla.

To nie są moje dane, a to wszystkie statystyki, które mam. Oczywiste jest, że rozkład nie jest normalny. Jedyne inne informacje, jakie posiadam, to to, że dane te reprezentują rządowe fundusze na mieszkańca dla różnych okręgów szkolnych.

Wiem wystarczająco dużo na temat statystyki, aby wiedzieć, że ten problem nie ma określonego rozwiązania, ale nie na tyle, aby wiedzieć, jak znaleźć właściwe domysły.

Czy rozkład logarytmiczny byłby odpowiedni? Jakich narzędzi mogę użyć do przeprowadzenia regresji (lub czy muszę to zrobić samodzielnie)?

Mark Eichenlaub
źródło
dodałem tag r, więc kod R jest podświetlony w moim komentarzu
mpiktas
Szczegółowy przykład tego samego pytania (i jego rozwiązania) można znaleźć w duplikacie wątku na stronie stats.stackexchange.com/questions/133129 .
whuber

Odpowiedzi:

17

Zastosowanie metody czysto statystycznej do wykonania tej pracy nie dostarczy absolutnie żadnych dodatkowych informacji na temat podziału wydatków szkolnych: wynik będzie jedynie odzwierciedlał arbitralny wybór algorytmu.

Potrzebujesz więcej danych .

Łatwo to osiągnąć: wykorzystuj dane z poprzednich lat, z porównywalnych dzielnic, cokolwiek. Na przykład federalne wydatki na 14866 okręgów szkolnych w 2008 r. Są dostępne na stronie spisu ludności . Pokazuje, że w całym kraju łączne dochody federalne na osobę (zarejestrowane) były w przybliżeniu logarytmicznie podzielone, ale podział według stanu wykazuje znaczne różnice ( np. Wydatki na kłody na Alasce mają ujemne przekrzywienie, podczas gdy wydatki na kłody w Kolorado mają silne dodatnie pochylenie) . Użyj tych danych, aby scharakteryzować prawdopodobną formę dystrybucji, a następnie dopasuj kwantyle do tej postaci.

Jeśli jesteś nawet blisko właściwej formy dystrybucyjnej, powinieneś być w stanie dokładnie odtworzyć kwantyle, dopasowując jeden lub co najwyżej dwa parametry. Najlepsza technika znajdowania dopasowania będzie zależeć od używanej formy dystrybucyjnej, ale - co ważniejsze - będzie zależeć od tego, do czego zamierzasz użyć wyników. Czy potrzebujesz oszacować średnią kwotę wydatków? Górne i dolne limity wydatków? Cokolwiek to jest, chcesz przyjąć pewną miarę dobroci dopasowania, która da ci najlepszą szansę na podejmowanie dobrych decyzji na podstawie wyników. Na przykład, jeśli twoje zainteresowanie koncentruje się na górnych 10% wszystkich wydatków, będziesz chciał dokładnie dopasować 95. percentyl i możesz nie dbać o dopasowanie 5. percentyla. Żadna wyrafinowana technika dopasowania nie sprawi, że te rozważania będą dla Ciebie.

Oczywiście nikt nie może w uzasadniony sposób zagwarantować, że ta oparta na danych, zorientowana na decyzje metoda będzie działała lepiej (lub gorzej) niż jakikolwiek przepis statystyczny, ale - w przeciwieństwie do podejścia czysto statystycznego - metoda ta ma podstawy oparte na rzeczywistości, koncentrując się na twoich potrzebach, nadając mu pewną wiarygodność i obronę przed krytyką.

Whuber
źródło
2
+1 Potrzebujesz więcej danych i tego, co zamierzasz wykorzystać, aby zasługiwać na dodatkowy nacisk.
vqv
2
To brzmi jak mądrość w twojej odpowiedzi. Będę musiał skonsultować więcej z ludźmi, którzy postawili mi problem dotyczący tego, czego chcą. Dziękuję za linki i porady.
Mark Eichenlaub
1
@ Mark Powodzenia!
whuber
23

Jak zauważył @whuber, metody statystyczne nie działają tutaj dokładnie. Musisz wywnioskować dystrybucję z innych źródeł. Kiedy znasz rozkład, masz ćwiczenie nieliniowego rozwiązywania równań. Oznacz przez funkcję kwantylową wybranego rozkładu prawdopodobieństwa za pomocą wektora parametru θ . Masz następujący nieliniowy układ równań:fθ

q0.05=f(0.05,θ)q0.5=f(0.5,θ)q0.95=f(0.95,θ)

gdzie to twoje kwantyle. Musisz rozwiązać ten system, aby znaleźć θ . Teraz praktycznie dla każdego rozkładu 3-parametrowego znajdziesz wartości parametrów spełniających to równanie. W przypadku rozkładów 2-parametrowych i 1-parametrowych ten system jest zawyżony, więc nie ma dokładnych rozwiązań. W takim przypadku możesz wyszukać zestaw parametrów, który minimalizuje rozbieżność:qθ

(q0.05f(0.05,θ))2+(q0.5f(0.5,θ))2+(q0.95f(0.95,θ))2

Tutaj wybrałem funkcję kwadratową, ale możesz wybrać, co chcesz. Zgodnie z komentarzami @whuber można przypisywać wagi, dzięki czemu ważniejsze kwantyle można dopasować dokładniej.

Dla czterech i więcej parametrów system jest niedookreślony, więc istnieje nieskończona liczba rozwiązań.

Oto przykładowy kod R ilustrujący to podejście. Dla celów demonstracyjnych generuję kwantyle z dystrybucji Singh-Maddala z pakietu VGAM . Ten rozkład ma 3 parametry i jest wykorzystywany w modelowaniu rozkładu dochodu.

 q <- qsinmad(c(0.05,0.5,0.95),2,1,4)
 plot(x<-seq(0,2,by=0.01), dsinmad(x, 2, 1, 4),type="l")
 points(p<-c(0.05, 0.5, 0.95), dsinmad(p, 2, 1, 4))

alternatywny tekst

Teraz utwórz funkcję, która ocenia nieliniowy układ równań:

 fn <- function(x,q) q-qsinmad(c(0.05, 0.5, 0.95), x[1], x[2], x[3])

Sprawdź, czy prawdziwe wartości spełniają równanie:

 > fn(c(2,1,4),q)
   [1] 0 0 0

Do rozwiązania układu równań nieliniowych używam funkcji nleqslvz pakietu nlqeslv .

 > sol <- nleqslv(c(2.4,1.5,4.3),fn,q=q)
 > sol$x       
  [1] 2.000000 1.000000 4.000001

Jak widzimy, otrzymujemy dokładne rozwiązanie. Spróbujmy teraz dopasować rozkład logarytmiczno-normalny do tych kwantyli. W tym celu użyjemy optimfunkcji.

 > ofn <- function(x,q)sum(abs(q-qlnorm(c(0.05,0.5,0.95),x[1],x[2]))^2)
 > osol <- optim(c(1,1),ofn)
 > osol$par
   [1] -0.905049  0.586334

Teraz wykreśl wynik

  plot(x,dlnorm(x,osol$par[1],osol$par[2]),type="l",col=2)
  lines(x,dsinmad(x,2,1,4))
  points(p,dsinmad(p,2,1,4))

alternatywny tekst

Z tego natychmiast widzimy, że funkcja kwadratowa nie jest tak dobra.

Mam nadzieję że to pomoże.

mpiktas
źródło
1
Świetny! Dzięki za cały wysiłek, który w to włożyłeś, mpiktas. Nie znam R, ale Twój kod jest wyjaśniony na tyle dobrze, że wciąż mogę łatwo powiedzieć, co robisz.
Mark Eichenlaub,
Wielkie dzięki za ten przykład. Myślę, że są 2 błędyofn <- function(x,q) sum(abs(q-qlnorm(c(0.05,0.5,0.95),x[1],x[2]))^2) . Proponuję, ofn <- function(x) sum(abs(q-qlnorm(c(0.05,0.5,0.95),x[1],x[2],x[3]))^2)ponieważ qnie jest wkładem ofni X[3]brakuje. Pozdrawiam
9

Wypróbuj pakiet rriskDistribution i - jeśli masz pewność co do lognormalnej rodziny dystrybucji - użyj polecenia

get.lnorm.par(p=c(0.05,0.5,0.95),q=c(8.135,11.259,23.611))

co powinno rozwiązać twój problem. Użyj fit.perczamiast tego, jeśli nie chcesz ograniczać się do jednego znanego pliku pdf.

Matthias Greiner
źródło
Super proste rozwiązanie!
luchonacho
6

Dla logarytmu normalnego stosunek 95 percentyla do mediany jest taki sam, jak stosunek mediany do 5 percentyla. To nawet nie jest prawie prawda, więc lognormal nie byłby dobrym pomysłem.

Masz wystarczająco dużo informacji, aby dopasować rozkład z trzema parametrami, i wyraźnie potrzebujesz rozkładu pochylenia. Dla analitycznej prostoty sugerowałbym przesunięty rozkład log-logistyczny jako jego funkcję kwantylową (tj. Odwrotność funkcji rozkładu skumulowanego) może być zapisana w dość prostej formie zamkniętej, więc powinieneś być w stanie uzyskać wyrażenia w postaci zamkniętej dla jego trzy parametry w kategoriach trzech kwantyli z odrobiną algebry (zostawię to jako ćwiczenie!). Ten rozkład jest wykorzystywany w analizie częstotliwości powodzi.

To jednak nie da ci żadnych wskazówek co do niepewności w szacunkach innych kwantyli. Nie wiem, czy jest to potrzebne, ale jako statystyk uważam, że powinienem móc to zapewnić, więc nie jestem zadowolony z tej odpowiedzi. Na pewno nie użyłbym tej metody lub prawdopodobnie żadnej metody do ekstrapolacji (znacznie) poza zakresem od 5 do 95 percentyli.

jeden przystanek
źródło
1
Dzięki za radę. Re: lognormal - Mógłbym ustalić stosunek percentyli do mediany, odejmując 7077 od wszystkiego, a następnie dodając go z powrotem na końcu. Jak zły byłby to pomysł?
Mark Eichenlaub
1
Dobrze, że dałoby to „przesunięty rozkład log-normalny”. Log-normal i log-logistic mają dość podobny kształt, z wyjątkiem cięższych ogonów tego drugiego, więc możesz wypróbować oba i porównać wyniki.
onestop
Porównaj jak? Przesunięty lognormal gwarantuje idealne dopasowanie do kwantyli. Prawie każda rodzina trzyparametrowa będzie idealnie pasować. Jak porównać dwa idealne dopasowania?
whuber
@ whuber Miałem na myśli porównanie uzyskanych prognoz dla percentyli odpowiadających innym wartościom
onestop
Coś mi brakuje: jakie inne wartości? PO stwierdza, że dostępne są tylko trzy percentyle, nic więcej.
whuber
2

Jedyne, co można wywnioskować z danych, to to, że rozkład jest niesymetryczny. Nie można nawet stwierdzić, czy te kwantyle pochodzą z dopasowanej dystrybucji, czy tylko z pliku ecdf.

Jeśli pochodzą one z dopasowanej dystrybucji, możesz wypróbować wszystkie dystrybucje, o których możesz pomyśleć i sprawdzić, czy pasują do siebie. Jeśli nie, nie ma prawie wystarczających informacji. Mógłbyś interpolować wielomian 2-go stopnia lub splajn 3-go stopnia dla funkcji kwantyli i użyć tego, lub wymyślić teorię dotyczącą rodziny rozkładów i dopasować kwantyle, ale wszelkie wnioski, jakie można by wyciągnąć za pomocą tych metod, byłyby głęboko podejrzane.

sesqu
źródło
1
Wielomiany i splajny prawdopodobnie nie będą prawidłowymi CDF.
whuber
Dobra obserwacja. W tym przypadku zwykły wielomian kwadratowy nie działa, ale istnieje nieskończenie wiele splajnów kwadratowych do wyboru (zdaniem Béziera), które nie powinny mieć tego samego problemu (chociaż niektóre mogą nadal wymagać przycinania domen). Podobnie powinno być możliwe znalezienie odpowiedniego monotonicznego splajnu sześciennego. Zdaję sobie sprawę z algorytmów splajnu, które gwarantują monotoniczność, ale nie jestem w stanie znaleźć takiego właśnie teraz, więc muszę zostawić sprawę „wybierz coś, co ci się podoba, działa jako cdf”.
sesqu
Możesz posunąć się tak daleko, aby dopasować monotoniczny splajn (lub cokolwiek innego) do logarytmów kwantyli, uzyskując w ten sposób coś rozsądnego w zakresie kwantyli. Ale to nie pomaga w dopasowaniu ogonów poza dwa skrajne kwantyle. Niechętnie należy pozostawić tak ważny aspekt dopasowania przypadkowej charakterystyce numerycznej procedury dopasowania.
whuber
2

Wykorzystanie kwantyli do oszacowania parametrów rozkładów a priori jest omówione w literaturze dotyczącej pomiaru czasu odpowiedzi człowieka jako „kwantylowe oszacowanie maksymalnego prawdopodobieństwa” (QMPE, choć pierwotnie błędnie nazwane „kwantylowym oszacowaniem maksymalnego prawdopodobieństwa”, QMLE), szczegółowo omówione przez Heathcote i koledzy . Możesz dopasować wiele różnych rozkładów a priori (wcześniej Gaussa, przesunięty Lognormal, Wald i Weibull), a następnie porównać prawdopodobieństwa logarytmu sumy wynikowych najlepszych dopasowań dla każdej dystrybucji, aby znaleźć smak dystrybucji, który wydaje się dawać najlepsze dopasowanie.

Mike Lawrence
źródło
2
Każdy rozkład trzyparametrowy gwarantuje idealne dopasowanie trzech kwantyli . Dlatego sensowne jest zastosowanie tego podejścia do dopasowania tylko jednego lub dwóch parametrów. Nie ma również sensu porównywanie dopasowania jednoparametrowego z dopasowaniem dwuparametrowym (z inną rodziną) na podstawie samego prawdopodobieństwa.
whuber
@ whuber, re: „Każdy rozkład trzyparametrowy gwarantuje idealne dopasowanie trzech kwantyli”. Nie zdawałem sobie z tego sprawy, więc dobrze wiedzieć! re: „Nie ma również sensu porównywanie dopasowania jednoparametrowego z dopasowaniem dwuparametrowym (z inną rodziną) na podstawie samego prawdopodobieństwa”. Ach tak, rzeczywiście; Nie wspomniałem, że należałoby zastosować pewną korektę złożoności (AIC, BIC, ...), jeśli porówna się dopasowania do smaków dystrybucji o różnej liczbie parametrów. Dzięki za zwrócenie na to uwagi.
Mike Lawrence
Troszkę przesadziłem, ponieważ myślałem o dwóch parametrach: skali i lokalizacji, a trzeci o szerokim zakresie kształtów. Mimo to większość rodzin trzyparametrowych ma wystarczającą elastyczność, aby zmieścić trzy percentyle, pod warunkiem że wszystkie są odrębne.
whuber
1

Możesz użyć informacji o percentylach, aby w pewien sposób zasymulować dane, a także użyć pakietu „logspline” pakietu R, aby oszacować rozkład nieparametryczny. Poniżej moja funkcja, która wykorzystuje taką metodę.

calc.dist.from.median.and.range <- function(m, r) 
{
    ## PURPOSE: Return a Log-Logspline Distribution given (m, r).
    ##          It may be necessary to call this function multiple times in order to get a satisfying distribution (from the plot). 
    ## ----------------------------------------------------------------------
    ## ARGUMENT:
    ##   m: Median
    ##   r: Range (a vector of two numbers)
    ## ----------------------------------------------------------------------
    ## RETURN: A log-logspline distribution object.
    ## ----------------------------------------------------------------------
    ## AUTHOR: Feiming Chen,  Date: 10 Feb 2016, 10:35

    if (m < r[1] || m > r[2] || r[1] > r[2]) stop("Misspecified Median and Range")

    mu <- log10(m)
    log.r <- log10(r)

    ## Simulate data that will have median of "mu" and range of "log.r"
    ## Distribution on the Left/Right: Simulate a Normal Distribution centered at "mu" and truncate the part above/below the "mu".
    ## May keep sample size intentionaly small so as to introduce uncertainty about the distribution. 
    d1 <- rnorm(n=200, mean=mu, sd=(mu - log.r[1])/3) # Assums 3*SD informs the bound
    d2 <- d1[d1 < mu]                   # Simulated Data to the Left of "mu"
    d3 <- rnorm(n=200, mean=mu, sd=(log.r[2] - mu)/3)
    d4 <- d3[d3 > mu]                   # Simulated Data to the Right of "mu"
    d5 <- c(d2, d4)                     # Combined Simulated Data for the unknown distribution

    require(logspline)
    ans <- logspline(x=d5)
    plot(ans)
    return(ans)
}
if (F) {                                # Unit Test 
    calc.dist.from.median.and.range(m=1e10, r=c(3.6e5, 3.1e12))
    my.dist <- calc.dist.from.median.and.range(m=1e7, r=c(7e2, 3e11))
    dlogspline(log10(c(7e2, 1e7, 3e11)), my.dist) # Density
    plogspline(log10(c(7e2, 1e7, 3e11)), my.dist) # Probability
    10^qlogspline(c(0.05, 0.5, 0.95), my.dist) # Quantiles 
    10^rlogspline(10, my.dist) # Random Sample 
}
Feiming Chen
źródło