Como estimar parâmetros para distribuição truncada Zipf a partir de uma amostra de dados?

Estou com um problema com o parâmetro de estimativa do Zipf. Minha situação é a seguinte:

Eu tenho um conjunto de amostras (medido a partir de um experimento que gera chamadas que devem seguir uma distribuição Zipf). Eu tenho que demonstrar que esse gerador realmente gera chamadas com distribuição zipf. Eu já li estas perguntas e respostas. Como calcular o coeficiente da lei de Zipf a partir de um conjunto de frequências principais? mas alcanço resultados ruins porque uso uma distribuição truncada. Por exemplo, se eu definir o valor "s" como "0,9" para o processo de geração, se tentar estimar o valor "s" conforme descrito nas perguntas e respostas relatadas, obtenho "s" igual a 0,2 ca. Eu acho que isso se deve ao fato de eu usar uma distribuição TRUNCATED (eu tenho que limitar o zipf com um ponto de truncamento, ele é truncado à direita).

Como posso estimar parâmetros com uma distribuição zipf truncada?

Atualização : 7 de abril de 2011 Esta resposta está ficando muito longa e cobre vários aspectos do problema em questão. No entanto, até agora resisti a dividi-lo em respostas separadas.

Adicionei na parte inferior uma discussão sobre o desempenho do de Pearson para este exemplo.$\chi^2$

Bruce M. Hill foi o autor, talvez, do artigo "seminal" sobre estimativa em um contexto semelhante ao Zipf. Ele escreveu vários trabalhos em meados da década de 1970 sobre o assunto. No entanto, o "estimador de Hill" (como agora é chamado) depende essencialmente das estatísticas de ordem máxima da amostra e, portanto, dependendo do tipo de truncamento presente, isso pode causar problemas.

O artigo principal é:

BM Hill, Uma abordagem geral simples à inferência sobre a cauda de uma distribuição , Ann. Estado. 1975.

Se seus dados forem inicialmente Zipf e, em seguida, truncados, uma boa correspondência entre a distribuição de graus e o gráfico Zipf poderá ser aproveitada.

Especificamente, a distribuição de graus é simplesmente a distribuição empírica do número de vezes que cada resposta inteira é vista,

Se traçarmos isso contra em um gráfico de log-log, obteremos uma tendência linear com uma inclinação correspondente ao coeficiente de escala.$i$

Por outro lado, se traçarmos o gráfico Zipf , onde classificaremos a amostra do maior para o menor e, em seguida, plotaremos os valores contra suas classificações, obteremos uma tendência linear diferente com uma inclinação diferente . No entanto, as pistas estão relacionadas.

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

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

$\hat{\beta}$

O @csgillespie deu um artigo recente, co-escrito por Mark Newman, em Michigan, sobre esse tópico. Ele parece publicar muitos artigos semelhantes sobre isso. Abaixo está outro, juntamente com algumas outras referências que podem ser interessantes. Newman às vezes não faz a coisa mais sensata estatisticamente, portanto, seja cauteloso.

MEJ Newman, Leis do poder, distribuições de Pareto e lei de Zipf , Contemporary Physics 46, 2005, pp. 323-351.

M. Mitzenmacher, Uma Breve História de Modelos Generativos para Lei de Potência e Distribuições Lognormal , Internet Math. vol. 1, n. 2, 2003, pp. 226-251.

K. Knight, Uma modificação simples do estimador de Hill com aplicações à robustez e redução de viés , 2010.

Adendo :

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

O gráfico resultante é

$i \leq 30$

Ainda assim, do ponto de vista prático, esse enredo deve ser relativamente atraente.

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

$\chi^2$

Também calcularemos uma segunda estatística formada pela primeira classificação das contagens em posições do tamanho 40, conforme mostrado na planilha de Maurizio (a última posição contém apenas a soma de vinte valores de resultados separados.

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

O papel

Clauset, A et al. , Distribuições da lei do poder em dados empíricos . 2009

contém uma descrição muito boa de como proceder para ajustar modelos de lei de energia. A página da web associada possui exemplos de código. Infelizmente, ele não fornece código para distribuições truncadas, mas pode fornecer um ponteiro.

Como um aparte, o artigo discute o fato de que muitos "conjuntos de dados de leis de poder" podem ser modelados igualmente bem (e, em alguns casos, melhor) com as distribuições normais ou exponenciais de Log!

Após a resposta detalhada do cardeal do usuário, realizei o teste do qui-quadrado na minha presumível distribuição zipf truncada. Os resultados do teste do qui-quadrado são relatados na tabela a seguir:

Onde StartInterval e EndInterval representam, por exemplo, o intervalo de chamadas e Observado é o número de chamadores gerando de 0 a 19 chamadas, etc. O teste do qui-quadrado é bom até o alcance das últimas colunas, elas aumentam cálculo, caso contrário, até esse ponto, o valor do qui-quadrado "parcial" era aceitável!

Com outros testes, o resultado é o mesmo, a última coluna (ou as duas últimas colunas) sempre aumenta o valor final e eu não sei por que e não sei se (e como) usar outro teste de validação.

PS: para ser completo, para calcular os valores esperados ( esperado ), sigo a sugestão do cardeal da seguinte maneira:

onde X_i 's são usados ​​para calcular :,x <- (1:n)^-S os P_i ' s para calcular p <- x / sum(x)e finalmente o E_i (número esperado de usuários para cada número de chamadas) é obtido porP_i * Total_Caller_Observed

e com Grau de liberdade = 13, a bondade do qui-quadrado rejeita sempre a hipótese de que o conjunto de amostras segue a distribuição Zipf porque a estatística de teste (64,14 neste caso) é maior do que a relatada nas tabelas do qui-quadrado, "demérito" para a última coluna. O resultado gráfico é relatado aqui:

embora o ponto de truncamento esteja definido como 500, o valor máximo obtido é 294. Penso que a "dispersão" final é a causa da falha do teste do qui-quadrado.

ATUALIZAR!!

Tento executar o teste do qui-quadrado em uma amostra de dados zipf presumível gerada com o código R relatado na resposta acima.

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

A plotagem associada é a seguinte:

Os resultados do teste do qui-quadrado são relatados na figura a seguir:

e a estatística do teste do qui-quadrado (44,57) é muito alta para a validação com o grau de liberdade escolhido. Também neste caso, a "dispersão" final dos dados é a causa do alto valor do qui-quadrado. Mas existe um procedimento para validar essa distribuição zipf (independentemente do meu gerador "errado", quero focar na amostra de dados R) ???