12

Tenho um grafo e calculei a distribuição dos graus e o grau como se segue:

dd <- degree_distribution(graph) 
 d <- degree(graph)

A partir disso, estimei a Power Law, para ver se a minha distribuição segue a "lei de Potência" :

    degree = 1:max(d)
    probability = dd[-1]

    # Exclui zeros, pois log de 0 é infinito!
    nonzero.position = which(probability != 0)
    probability = probability[nonzero.position]
    degree = degree[nonzero.position]

    reg = lm(log(probability) ~ log(degree))
    cozf = coef(reg)
    #Estima a power law com base nos valores 
    power.law.fit = function(x) exp(cozf[[1]] + cozf[[2]] * log(x))

A partir disso, plotei os pontos e a power law usando o ggplot2.
Resultando na imagem a seguir:

df <- data.frame(x = degree, y = probability)
  print(
      ggplot(df, aes(x,y,colour="Distribuição"))+
        geom_point(shape = 4) +
        stat_function(fun = power.law.fit, geom = "line", aes(colour="Power Law"))+

        labs(title = "Grafo", subtitle = "Distribuição dos Graus",
             x="K", y="P(k)", colour="Legenda")+
        scale_color_brewer(palette="Dark2")
  )

inserir a descrição da imagem aqui

Como é possível observar, a minha distribuição não segue a Power Law! Gostaria de estimar a distribuição de Poisson e plotar no mesmo gráfico.
Mesmo não tendo certeza que a minha distribuição não segue (ou segue) a Poisson, gostaria de plotar junto com a Power Law. Não tenho a mínima ideia de como estimar essa distribuição (Poisson) a partir dos dados, e calcular o grau médio.

Alguém poderia me ajudar? Grato.

  • O grafo usado para calcular a distribuição e o grau é muito grande (700 mil vértices), por isso não coloquei os dados dos grafos. A explicação da resposta pode ser baseada em um grafo qualquer.
5
  • Filipe, apesar de a pergunta estar bem formulada, ela parece fugir do escopo do site infelizmente, pois trata mais de estatística do que de programação. Commented 31/08/2017 às 8:30
  • 7
    Não acredito que fuja do escopo, estatística faz parte da programação, e a própria tag r diz R é um ambiente e linguagem de programação de código aberto para computação estatística...
    – Raizant
    Commented 31/08/2017 às 13:54
  • Também concordo @knautiluz. Espero que alguém me ajude kk. Agradecido.
    – Fillipe
    Commented 31/08/2017 às 16:02
  • @Knautiluz estatística não faz parte da programação. Perguntas de estatística estão fora do escopo do site. Veja mais aqui pt.meta.stackoverflow.com/questions/4266/… Commented 1/09/2017 às 17:15
  • 1
    Ele perguntou usando linguagem r, que é uma linguagem de programação voltada para estatística, ainda que a pergunta seja mais voltada para estatística. Na maioria das faculdades de programação nós temos que aprender estatística para poder lidar com certos casos envolvendo programação. Eu entendi o que você disse, mas mesmo assim não acredito que fuja do escopo. Agora sobre o Cross Validated sim, realmente seria melhor respondida lá, acho isso poderia ter sido comentado de primeira.
    – Raizant
    Commented 1/09/2017 às 18:24

1 Resposta 1

10
+50

Existem várias formas de estimar os parâmetros de uma distribuição (máxima verossimilhança, métodos de momentos, Bayes). Entrar nesse âmbito fugiria ao escopo do site, pois é uma pergunta de estatística --- isso seria melhor respondido no cross-validated.

Dito isto, você pode estimar por máxima verossimilhança usando o pacote MASS. Suponha que os seus dados estejam na variável x, vou simular uns dados para o exemplo:

rm(list = ls())
set.seed(10)
x <- rpois(n = 100, lambda = 10)

A função fitdistr do pacote MASS faz o ajuste por máxima verossimilhança (você também pode calcular facilmente na mão, derivando a função de log-verossimilhança e otimizando a função):

library(MASS)
lambda <- fitdistr(x, "poisson")
lambda
  10.2700000 
 ( 0.3204684)

Veja que o valor estimado ficou bem próximo do valor real.

Você também poderia ter estimado com a função de modelos lineares generalizados do R:

glm(x ~ 1, family = poisson(link = identity))

O que vai te dar o mesmo resultado.

Entretanto, para o caso da distribuição de Poisson você nem precisa se dar ao trabalho de fazer isso. A estimativa de máxima verossimilhança do parâmetro lambda da distribuição de poisson é simplesmente a média:

mean(x)
[1] 10.27

Uma vez que você tem o valor de lambda, você pode calcular a densidade da distribuição de poisson usando a função dpois() e colocar esses valores no seu gráfico. Um exemplo abaixo:

hist(x, freq = FALSE, col = "lightblue")
seq <- seq.int(0, max(x))
lines(seq, dpois(seq, lambda = lambda$estimate), col = "red")

inserir a descrição da imagem aqui

Para fazer o gráfico com ggplot2:

library(ggplot2)
df  <- data.frame(x = x)
df2 <- data.frame(seq, dens = dpois(seq, lambda = mean(x)))
ggplot(df, aes(x = x)) + 
  geom_histogram(aes(y = ..density..), binwidth = 1, col = "black", fill = "lightblue") +
  geom_bar(aes(x = seq, y = dens), data = df2, col = "red", lwd = 1.2, width = 0.0001,  stat = "identity")

inserir a descrição da imagem aqui

Note que eu mudei para barras ao invés de linha, pois a distribuição de poisson é discreta. Mas se você quiser colocar linhas basta mudar geom_bar() para geom_line(). Para fazer essa alteração no gráfico base, adicione o parâmetro type = "h" na função lines().

7
  • como eu plotaria usando o ggplot2?
    – Fillipe
    Commented 2/09/2017 às 18:31
  • @Fillipe coloquei o exemplo do ggplot2 Commented 3/09/2017 às 2:53
  • fiz o seguinte: seq <- seq.int(0, max(dgDist)) df2 <- data.frame(seq, dens = dpois(seq, lambda = mean(dgDist))) dgDist é a minha distribuição dos graus. Usei: geom_line(aes(x = seq, y = dens, colour="Poisson"), data = df2) mas não plota a Poisson, só plota os dados.
    – Fillipe
    Commented 3/09/2017 às 15:26
  • Percebi que ele não está criando uma sequência quando uso seq <- seq(0, max(dgDist)) (pois o maior valor é 0.75). Usei seq <- seq(0, max(dgDist),0.05) e ele plotou um gráfico assim: prntscr.com/ggkfvh
    – Fillipe
    Commented 3/09/2017 às 16:01
  • @Fillipe não tem como o número maximo ser 0.75, pois a variável é discreta, você está usando a variável errada. Commented 3/09/2017 às 17:01

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .