估计泊松分布

Question

我有一个图表，我计算了度数分布和度数如下：

library(igraph) # for these two functions

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

由此，我估计了幂律，看看我的分布是否遵循“幂律”：

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

nonzero.position = which(probability != 0)
probability = probability[nonzero.position]
degree = degree[nonzero.position]

reg = lm(log(probability) ~ log(degree))
cozf = coef(reg)

power.law.fit = function(x) exp(cozf[[1]] + cozf[[2]] * log(x))

由此，我使用ggplot2 绘制了点和幂律。 结果如下图：

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

如您所见，我的分布不遵循幂律！我想估计泊松分布并绘制在同一张图上。 尽管我不确定我的分布不遵循（或遵循）泊松，但我想将它与幂律一起绘制。我不知道如何从数据中估计这个分布（泊松），并计算平均度数。

谁能帮帮我？

用于计算分布和度数的图非常大（70万个顶点），所以我没有放图的数据。答案的解释可以基于任何图表。

Answer 1

来自?dpois：

泊松分布有密度

p(x) = λ^x exp(-λ)/x!

对于 x = 0, 1, 2, ... 。均值和方差为E(X) = Var(X) = λ。

所以我将使用秘密 lambda 生成一些虚拟数据：

mysecret <- ####

x <- data.frame(xes = rpois(50, mysecret))

> x$xes
 [1] 0 2 2 1 1 4 1 1 0 2 2 2 1 0 0 1 2 3 2 4 2 1 0 3 2 1 3 1 2 1 5 0 2 3 2 1 0 1 2 3 0 1 2 2 0 3 2 2 2 3


> mean(x$xes)
[1] 1.66
> var(x$xes)
[1] 1.371837

所以我的秘密 lambda 有两个很好的猜测是 1.66 和 1.37。让我们试试吧：

library(ggplot2)
ggplot(x, aes(xes)) + 
  geom_histogram(aes(y = ..density.., color = "Raw data"), 
                 fill = "white", binwidth = 1, center = 0, size = 1.5) +
  stat_summary(fun.y = dpois, aes(x = xes, y = xes, color = "Density based on E(X)"), 
               fun.args = list(lambda = 1.66), geom = "line", size = 1.5) +
  stat_summary(fun.y = dpois, aes(x = xes, y = xes, color = "Density based on Var(X)"), 
               fun.args = list(lambda = 1.37), geom = "line", size = 1.5)

他们俩都很好。您不能真正使用内置的stat_function 或geom_density 来生成这些，因为泊松分布仅针对整数定义。直方图和汇总函数运行良好，因为它们仅在数据点本身进行估计，而不是插值。

如果您想了解更多细节，可以使用MASS 包：

MASS::fitdistr(x$xes, dpois, start = list(lambda = 1))

    lambda  
  1.6601563 
 (0.1822258)

所以让我们尝试从它构建：

library(dplyr)
df <- data_frame(xes = seq.int(max(x$xes)+1)-1,
                 dens.m = dpois(xes, 1.66),
                 dens.u = dpois(xes, 1.66+0.18),
                 dens.l = dpois(xes, 1.66-0.18))

> df
# A tibble: 6 x 4
    xes     dens.m     dens.u     dens.l
  <dbl>      <dbl>      <dbl>      <dbl>
1     0 0.19013898 0.15881743 0.22763769
2     1 0.31563071 0.29222406 0.33690378
3     2 0.26197349 0.26884614 0.24930880
4     3 0.14495866 0.16489230 0.12299234
5     4 0.06015785 0.07585046 0.04550717
6     5 0.01997240 0.02791297 0.01347012

ggplot(x, aes(xes)) + 
  geom_histogram(aes(y = ..density..), color = "black",
                 fill = "white", binwidth = 1, center = 0, size = 1.5) +
  geom_ribbon(data = df, aes(xes, ymin = dens.l, ymax = dens.u), fill = "grey50", alpha = 0.5) +
  geom_line(data = df, aes(xes, dens.m, color = "Based on E(X)\n+/-1 SD of lambda"), size = 1.5)

基于这两种方法和视觉解释，您应该可以轻松地说 λ = 1.66+/-0.18。

作为参考，我的秘密初始值为 1.5。