Jak oszacować parametry dla obciętej dystrybucji Zipf na podstawie próbki danych?

Mam problem z parametrem oszacowania dla Zipf. Moja sytuacja jest następująca:

Mam zestaw próbek (mierzony na podstawie eksperymentu, który generuje połączenia, które powinny być zgodne z rozkładem Zipf). Muszę wykazać, że ten generator naprawdę generuje połączenia z dystrybucją zipf. Przeczytałem już to pytanie. Jak obliczyć współczynnik prawa Zipfa z zestawu najwyższych częstotliwości? ale osiągam złe wyniki, ponieważ używam skróconej dystrybucji. Na przykład, jeśli ustawię wartość „s” na „0,9” w procesie generowania, to jeśli spróbuję oszacować wartość „s”, jak napisano w zgłoszonych pytaniach i odpowiedziach, otrzymam „s” równe 0,2 ca. Myślę, że wynika to z faktu, że używam dystrybucji TRUNCATED (muszę ograniczyć zipf z punktem obcięcia, jest on skrócony w prawo).

Jak mogę oszacować parametry przy obciętym rozkładzie zipf?

distributions estimation pareto-distribution zipf Maurizio
źródło

dla jasności, co dokładnie obciąłeś? Rozkład wartości czy sam wykres Zipf? Czy znasz punkt obcięcia? Czy obcięcie jest artefaktem danych czy artefaktem przetwarzania danych (np. Jakąś decyzję podjętą przez ciebie lub eksperymentatora)? Pomocne będą wszelkie dodatkowe szczegóły.

kardynał

@kardynał. (część 1/2) Dzięki kardynałowi. Podam więcej szczegółów: Mam generator VoIP, który generuje połączenia po Zipf (i innej dystrybucji) dla głośności na rozmówcę. Muszę sprawdzić, czy ten generator naprawdę podąża za tymi dystrybucjami. W przypadku Zipf Distribution muszę zdefiniować punkt obcięcia (stąd jest znany i odnosi się do rozkładu wartości), który jest maksymalną liczbą wygenerowanych połączeń przez użytkownika i parametru skali. W szczególności w moim przypadku wartość ta wynosi 500, co oznacza, że jeden użytkownik może wygenerować maksymalnie 500 połączeń.

Maurizio

(część 2/2) Innym ustawianym parametrem jest parametr skali dla Zipf, który określa rozkład rozkładu (ta wartość w moim przypadku wynosi 0,9). Mam wszystkie parametry (wielkość próbki, częstotliwość na użytkownika itp.), Ale muszę sprawdzić, czy mój zestaw danych jest zgodny z dystrybucją zipf.

Maurizio

więc najwyraźniej renormalizujesz rozkład o

, ponieważ dla tego, co uważałbym za „obcięty Zipf”, parametr skalowania 0,9 byłby niemożliwy. Jeśli możesz wygenerować wiele tych danych i „tylko” masz 500 możliwych wyników, dlaczego nie skorzystać z testu dobroci dopasowania chi-kwadrat? Ponieważ twoja dystrybucja ma długi ogon, możesz potrzebować dość dużej próby. Ale to byłby jeden sposób. Inną szybką i brudną metodą byłoby sprawdzenie, czy otrzymujesz właściwy rozkład empiryczny dla małych wartości liczby połączeń.

\sum_{i = 1}^{500} i^{- 0.9}

$\sum_{i=1}^{500} i^{-0.9}$

kardynał

Odpowiedzi:

Aktualizacja : 7 kwietnia 2011 r. Ta odpowiedź jest dość długa i obejmuje wiele aspektów aktualnego problemu. Jednak jak dotąd opierałem się, dzieląc go na osobne odpowiedzi.

Na samym dole dodałem dyskusję na temat wydajności Pearsona dla tego przykładu. $\chi^2$

Być może Bruce M. Hill jest autorem „przełomowego” artykułu na temat szacunków w kontekście podobnym do Zipf. W połowie lat 70. napisał na ten temat kilka artykułów. Jednak „estymator Hill” (jak się teraz nazywa) zasadniczo opiera się na statystykach maksymalnego rzędu próbki, a zatem, w zależności od rodzaju obcięcia, może sprawić ci kłopotów.

Główny artykuł to:

BM Hill, Proste ogólne podejście do wnioskowania na temat ogona dystrybucji , Ann. Stat. , 1975.

Jeśli twoje dane naprawdę są początkowo Zipf, a następnie są obcinane, to dobra korespondencja między rozkładem stopni a działką Zipf może zostać wykorzystana na Twoją korzyść.

W szczególności rozkład stopni jest po prostu rozkładem empirycznym liczby wyświetleń każdej odpowiedzi całkowitej,

d_{i} = \frac{# {j : X_{j} = i}}{n} .

$d_i = \frac{\#\{j: X_j = i\}}{n} .$

Jeśli narysujemy to względem na wykresie log-log, otrzymamy trend liniowy o nachyleniu odpowiadającym współczynnikowi skalowania. $i$

Z drugiej strony, jeśli wykreślimy wykres Zipf , w którym sortujemy próbkę od największej do najmniejszej, a następnie wykreślamy wartości względem ich rang, otrzymujemy inny trend liniowy z innym nachyleniem. Jednak stoki są powiązane.

$\alpha$ $-\alpha$ $-1/(\alpha-1)$ $\alpha = 2$ $n = 10^6$ $-2$ $-1/(2-1) = -1$

Wykresy rozkładu stopni (po lewej) i Zipf (po prawej) dla próbki iid z rozkładu Zipf.

$\tau$ $\tau$ $\alpha$

$\hat{\beta}$

\hat{α} = 1 - \frac{1}{\hat{β}} .

$\hat{\alpha} = 1 - \frac{1}{\hat{\beta}} .$

@csgillespie opublikował jeden z ostatnich artykułów współautora Marka Newmana z Michigan na ten temat. Wydaje się, że publikuje wiele podobnych artykułów na ten temat. Poniżej znajduje się kolejna wraz z kilkoma innymi referencjami, które mogą być interesujące. Newman czasami nie robi statystycznie najbardziej sensownej rzeczy, więc bądź ostrożny.

MEJ Newman, Prawa potęgi, rozkłady Pareto i prawo Zipfa , Contemporary Physics 46, 2005, s. 323–351.

M. Mitzenmacher, Krótka historia modeli generatywnych dla prawa mocy i rozkładów logarytmicznych , matematyka internetowa. , vol. 1, nr 2, 2003, s. 226–251.

K. Knight, Prosta modyfikacja estymatora Hill'a z aplikacjami do odporności i redukcji uprzedzeń , 2010.

Dodatek :

$R$ $10^5$

> x <- (1:500)^(-0.9)
> p <- x / sum(x)
> y <- sample(length(p), size=100000, repl=TRUE, prob=p)
> tab <- table(y)
> plot( 1:500, tab/sum(tab), log="xy", pch=20, 
        main="'Truncated' Zipf simulation (truncated at i=500)",
        xlab="Response", ylab="Probability" )
> lines(p, col="red", lwd=2)

Powstały wykres to

„Skrócony” wykres Zipf (obcięty przy i = 500)

$i \leq 30$

Jednak z praktycznego punktu widzenia taka fabuła powinna być względnie atrakcyjna.

$\alpha = 2$ $n = 300\,000$ $x_{\mathrm{max}} = 500$

$\chi^2$

X^{2} = \sum_{i = 1}^{500} \frac{(O_{i} - E_{i})^{2}}{E_{i}}

$X^2 = \sum_{i=1}^{500} \frac{(O_i - E_i)^2}{E_i}$

O_{i}

$O_i$

i

$i$

E_{i} = n p_{i} = n i^{- α} / \sum_{j = 1}^{500} j^{- α}

$E_i = n p_i = n i^{-\alpha} / \sum_{j=1}^{500} j^{-\alpha}$

Obliczymy również drugą statystykę utworzoną przez pierwsze binowanie liczb w pojemnikach o rozmiarze 40, jak pokazano w arkuszu kalkulacyjnym Maurizio (ostatni bin zawiera tylko sumę dwudziestu oddzielnych wartości wyników.

$n$ $p$

$p$

wprowadź opis zdjęcia tutaj

$R$

# Chi-square testing of the truncated Zipf.

a <- 2
n <- 300000
xmax <- 500

nreps <- 5000

zipf.chisq.test <- function(n, a=0.9, xmax=500, bin.size = 40)
{
  # Make the probability vector
  x <- (1:xmax)^(-a)
  p <- x / sum(x)

  # Do the sampling
  y <- sample(length(p), size=n, repl=TRUE, prob=p)

  # Use tabulate, NOT table!
  tab <- tabulate(y,xmax)

  # unbinned chi-square stat and p-value
  discrepancy <- (tab-n*p)^2/(n*p)
  chi.stat <- sum(discrepancy)
  p.val    <- pchisq(chi.stat, df=xmax-1, lower.tail = FALSE)

  # binned chi-square stat and p-value
  bins <- seq(bin.size,xmax,by=bin.size)
  if( bins[length(bins)] != xmax )
    bins <- c(bins, xmax)

  tab.bin  <- cumsum(tab)[bins]
  tab.bin <- c(tab.bin[1], diff(tab.bin))

  prob.bin <- cumsum(p)[bins] 
  prob.bin <- c(prob.bin[1], diff(prob.bin))

  disc.bin <- (tab.bin - n*prob.bin)^2/(n * prob.bin)
  chi.stat.bin <- sum(disc.bin)
  p.val.bin <- pchisq(chi.stat.bin, df=length(tab.bin)-1, lower.tail = FALSE)

  # Return the binned and unbineed p-values
  c(p.val, p.val.bin, chi.stat, chi.stat.bin)
}

set.seed( .Random.seed[2] )

all <- replicate(nreps, zipf.chisq.test(n, a, xmax))

par(mfrow=c(2,1))
hist( all[1,], breaks=20, col="darkgrey", border="white",
      main="Histogram of unbinned chi-square p-values", xlab="p-value")
hist( all[2,], breaks=20, col="darkgrey", border="white",
      main="Histogram of binned chi-square p-values", xlab="p-value" )

type.one.error <- rowMeans( all[1:2,] < 0.05 )

kardynał
źródło

+1, jak zwykle świetna odpowiedź. Powinieneś wyznaczyć siebie na moderatora, pozostała jeszcze 1 godzina :)

mpiktas

@mpiktas, dzięki za komplementy i zachęty. Nie jestem pewien, czy mógłbym uzasadnić mianowanie się i tak już bardzo silną listą kandydatów, którzy jednakowo uczestniczyli w większym i dłuższym okresie niż ja.

kardynał

@cardinal, oto kilka linków do alternatywa dla estymatora Hilla: Oryginalny artykuł przez Paulauskas i następują-ups by Vaiciulis i Gadeikis i Paulauskas . Ten estymator miał prawdopodobnie lepsze właściwości niż oryginalne Hill's.

mpiktas

@mpiktas, dzięki za linki. Istnieje całkiem sporo „nowych i ulepszonych” wersji estymatora Hill. Główną wadą oryginalnego podejścia jest to, że wymaga ono wyboru „punktu odcięcia” w miejscu, gdzie zatrzymać uśrednianie. Myślę, że w większości zostało to zrobione przez „spoglądanie w oczy”, co otwiera nas na zarzut podmiotowości. Jedna z książek Resnicka na temat dystrybucji o długich ogonach omawia to bardziej szczegółowo, o ile pamiętam. Myślę, że to jego nowsza wersja.

kardynał

@cardinal, dziękuję bardzo, jesteś bardzo miły i bardzo szczegółowy! Twój przykład w R był dla mnie bardzo przydatny, ale jak mogę w tym przypadku przeprowadzić formalny test chi-kwadrat? (Użyłem testu chi-kwadrat z innymi rozkładami, takimi jak jednolity, wykładniczy, normalny, ale mam wiele wątpliwości co do zipf .. Przepraszam, ale to moje pierwsze podejście do tych tematów). Pytanie do modetatorów: czy muszę napisać kolejne pytania i odpowiedzi, takie jak „jak wykonać test chi-kwadrat dla okrojonej dystrybucji zipf?” lub kontynuować w pytaniach i odpowiedziach, może aktualizując tagi i tytuł?

Maurizio

Papier

Clauset, A i in. , Power-law Distribution in Empirical Data . 2009

zawiera bardzo dobry opis sposobu dopasowania modeli prawa mocy. Powiązana strona internetowa zawiera próbki kodu. Niestety nie podaje kodu dla skróconych dystrybucji, ale może dać ci wskaźnik.

Nawiasem mówiąc, w artykule omówiono fakt, że wiele „zestawów danych dotyczących prawa mocy” można modelować równie dobrze (a w niektórych przypadkach lepiej) za pomocą rozkładu normalnego lub wykładniczego Log!

csgillespie
źródło

Niestety, ten artykuł nie mówi nic o skróconej dystrybucji. Znalazłem kilka pakietów w R, które radzą sobie z parametrem estymacji Zipf w prosty sposób (zipfR, VGAM), ale skrócona dystrybucja wymaga „specjalnego traktowania”. Czy w ostatnim zdaniu miałeś na myśli, że możliwe jest modelowanie zbioru danych dotyczących mocy z np. Rozkładem wykładniczym, a następnie zastosowanie procesu parametrów estymacji dla „obciętego” rozkładu wykładniczego? Jestem bardzo nowy w tym temacie!

Maurizio

W artykule autorzy ponownie analizują różne zestawy danych, w których wprowadzono prawo mocy. Autorzy zwracają uwagę, że w wielu przypadkach model prawa mocy nie jest tak świetny, a alternatywna dystrybucja byłaby lepsza.

csgillespie

Po szczegółowej odpowiedzi kardynała użytkownika wykonałem test chi-kwadrat na moim przypuszczalnie obciętym rozkładzie zipf. Wyniki testu chi-kwadrat podano w poniższej tabeli:

wprowadź opis zdjęcia tutaj

Tam, gdzie StartInterval i EndInterval reprezentują na przykład zakres połączeń, a Obserwowana to liczba dzwoniących generujących od 0 do 19 połączeń itd. Test chi-kwadrat jest dobry do osiągnięcia ostatnich kolumn, zwiększają końcową obliczenia, w przeciwnym razie do tego momentu akceptowalna była „częściowa” wartość chi-kwadrat!

W przypadku innych testów wynik jest taki sam, ostatnia kolumna (lub ostatnie 2 kolumny) zawsze zwiększa wartość końcową i nie wiem, dlaczego i nie wiem, czy (i jak) użyć innego testu sprawdzania poprawności.

PS: dla kompletności, aby obliczyć oczekiwane wartości ( oczekiwane ), postępuję zgodnie z sugestią kardynała w ten sposób:

wprowadź opis zdjęcia tutaj

gdzie x_i „s są wykorzystywane do obliczenia: x <- (1:n)^-SThe P_i ” s do obliczania p <- x / sum(x)i wreszcie E_i (Oczekiwany nr użytkowników dla każdego nr połączeń) otrzymuje sięP_i * Total_Caller_Observed

a przy stopniu swobody = 13 dobroć chi-kwadrat zawsze odrzuca hipotezę, że zestaw próbek jest zgodny z rozkładem Zipf, ponieważ statystyki testowe (w tym przypadku 64,14) są większe niż te podane w tabelach chi-kwadrat „demerit” dla ostatniej kolumny. Wynik graficzny przedstawiono tutaj: wprowadź opis zdjęcia tutaj

chociaż punkt obcięcia jest ustawiony na 500, maksymalna uzyskana wartość to 294. Myślę, że ostateczna „dyspersja” jest przyczyną niepowodzenia testu chi-kwadrat.

AKTUALIZACJA!!

Próbuję wykonać test chi-kwadrat na przypuszczalnej próbce danych zipf wygenerowanej za pomocą kodu R podanego w odpowiedzi powyżej.

> x <- (1:500)^(-2)
> p <- x / sum(x)
> y <- sample(length(p), size=300000, repl=TRUE, prob=p)
> tab <- table(y)
> length(tab)
[1] 438
> plot( 1:438, tab/sum(tab), log="xy", pch=20, 
        main="'Truncated' Zipf simulation (truncated at i=500)",
        xlab="Response", ylab="Probability" )
> lines(p, col="red", lwd=2)

Powiązana fabuła jest następująca: wprowadź opis zdjęcia tutaj

Wyniki testu chi-kwadrat przedstawiono na poniższym rysunku: wprowadź opis zdjęcia tutaj

a statystyka testu chi-kwadrat (44,57) jest zbyt wysoka, aby można było przeprowadzić walidację przy wybranym stopniu wolności. Również w tym przypadku ostateczne „rozproszenie” danych jest przyczyną wysokiej wartości chi-kwadrat. Ale istnieje procedura sprawdzania poprawności dystrybucji zipf (niezależnie od mojego „złego” generatora, chcę skupić się na próbce danych R)?

Maurizio
źródło

@ Maurizio, z jakiegoś powodu brakowało mi tego postu do tej pory. Czy w każdym razie możesz go edytować i dodać fabułę podobną do ostatniej w moim poście, ale używając obserwowanych danych? To może pomóc zdiagnozować problem. Wydaje mi się, że widziałem inne pytanie, w którym miałeś problem z uzyskaniem jednolitego rozkładu, więc może to również dotyczy tych analiz. (?) Pozdrowienia.

kardynał

@cardinal, zaktualizowałem wyniki! Co myślisz? Pytanie o rozkład równomierny to kolejna rzecz, którą muszę sprecyzować lepiej i zrobię to dzisiaj lub jutro;)

Maurizio,

S = 0.9

$S = 0.9$

p = P (X_{i} = 500) \approx 4.05 \times 10^{- 4}

$p = \mathbb{P}(X_i = 500) \approx 4.05 \times 10^{-4}$

n = 8454

$n = 8454$

8454 \cdot 4.05 \cdot 10^{- 4} \approx 3.43

$8454 \cdot 4.05 \cdot 10^{-4} \approx 3.43$

1 - (1 - 0.000405)^{8454} \approx 0.9675

$1 - (1 - 0.000405)^{8454} \approx 0.9675$ . Zwróć uwagę, jak blisko to pasuje do powyższej symulacji.

kardynał

@ cardinal, myślę też, że jest coś „nie tak” w procedurze generowania (moim celem jest sprawdzenie, czy ten generator naprawdę podąża za dystrybucją Zipf). W tych dniach muszę rozmawiać z projektantami projektu.

Maurizio