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?

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,

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$

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

$\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

$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$

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$

$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 )

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!

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:

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:

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:

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:

Wyniki testu chi-kwadrat przedstawiono na poniższym rysunku:

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)?