利用Log Binning绘图拟合参数

发表于 2022-04-05 更新于 2022-04-12 分类于绘图本文字数： 8.1k 阅读时长 ≈ 7 分钟

绘制度分布是分析网络属性的一个组成部分。该过程从获得\(N_{k}\)开始，即度数为\(k\)的节点数。这可以通过直接测量或模型来提供。从\(N_{k}\)我们计算出\(p_{k}=N_{k}/N\)。问题是，如何绘制\(p_{k}\)以最好地提取其属性。

使用log-log图

在无标度网络中，具有一或两条链路的众多节点与少数节点共存，其中少数节点为具有数千甚至数百万链路的节点。使用线性 k 轴压缩无数小k区域中的节点，使它们不可见。类似地，由于\(k=1\)和大\(k\)的\(p_{k}\)可能存在数量级差异，如果我们在线性垂直轴上绘制\(p_{k}\)，大\(k\)的值将显示为零（图 4.22a）。对数图的使用避免了这些问题。我们可以使用10次方的对数轴（图 4.22b），或者我们可以绘制\(\log k\)函数的\(\log k\)。请注意，\(p_{k}=0\)或\(k=0\)的点不会在\(\log - \log\)图上显示，因为\(\log0=-\infty\)。

避免Linear Binning

最有缺陷的方法（但在文献中经常出现）是在对数图上简单地绘制\(p_{k}=N_{k}/N\)（图 4.22b）。这称为线性分箱（Linear Binning），因为每个bin具有相同的大小\(\Delta k=1\)。对于无标度网络，linear binning会在大\(k\)处产生显而易见的平台，由形成水平线的大量数据点组成（图 4.22b）。这个平台有一个简单的解释：通常我们只有一个高度节点的样本，因此在高 \(k\)区域中，我们要么有\(N_{k}=0\)（没有具有\(k\)度的节点），要么有\(N_{k}=1\)（具有\(k\)度的单个节点）。因此，Linear Binning将提供\(p_{k}=0\)（未在对数图上显示）或\(p_{k}=1/N\)（适用于所有hubs），在\(p_{k}=1/N\)处生成一个平台。这个平台会影响我们估计度指数\(\gamma\)的能力。例如，如果我们尝试使用linear binning对图4.22b中所示的数据拟合幂律，则获得的\(\gamma\)与实际值\(\gamma =2.5\)完全不同。原因是在linear binning下，我们在小\(k\)的bin中有大量节点，这使我们能够自信地在这种情况下拟合\(p_{k}\)。在大\(k\)的bin中，我们的节点太少，无法对\(p_{k}\)进行适当的统计估计。相反，新出现的平台会使得拟合参数偏离。然而，正是这种高\(k\)状态在确定\(\gamma\)中起关键作用。增加bin大小不会解决这个问题。因此，建议避免对肥尾分布进行Linear binning。

使用Logarithmic Binning

Logarithmic binning纠正了linear binning的非均匀采样。对于log-binning，我们让bin大小随程度增加，确保每个bin具有相当数量的节点。例如，我们可以选择bin大小为2的倍数，这样第一个bin的大小为\(b_{0}=1\)，包含所有\(k=1\)的节点；第二个大小为\(b_{1}=2\)，包含度数\(k=2,3\)的节点；第三个bin的大小为\(b_{2}=4\)，包含度数\(k=4,5,6,7\)的节点。通过归纳，第\(n\)个bin的大小为\(2^{n-1}\)，包含度数为\(k=2^{n-1},2^{n-1}+1,...,2^{n}-1\)的节点。请注意，bin大小可以随任意增量增加，\(b_{n}=c^{n}\)，其中\(c>1\)。度分布由\(p_{\left \langle k \right \rangle}=N_{n}/\left(Nb_{n} \right)\)给出，其中\(N_{n}\)是在大小为\(b_{n}\)的\(bin_{n}\)中找到的节点数，\(\left \langle k_{n} \right \rangle\)是bin\(b_{n}\)中节点的平均度数。图4.22c显示了logarithmic binning的\(p_{k}\)。请注意，现在扩展到高\(k\)平台，其本来在linear binning下不可见。因此，logarithmic binning也可以从稀有的高度节点中提取有用信息。由于上述操作相当于把每个bin中的度的\(p_{k}\)进行平均，所以最终在高\(k\)的bin中有些\(p_{k}\)是0，所以平均之后的值要小于\(p_{k}=1/N\)，这是要值得注意的。

使用累积分布（Cumulative Distribution）

从\(p_{k}\)的尾部提取信息的另一种方法是绘制互补累积分布， \[ \begin{equation} P_{k}=\sum^{\infty}_{q=k+1}p_{q}, \end{equation} \] 这再次增强了高k区域的统计显著性。如果\(p_{k}\)遵循幂律\(p_{k}=k^{-\gamma}\)，则累积分布缩放为 \[ \begin{equation} p_{k}\sim k^{-\gamma+1}. \end{equation} \] 累积分布再次消除了linear binning观察到的平台,并扩展了区域（图 4.22d），从而可以更准确地估计度指数。

注意

当横坐标并非为离散的变量时，需要先把连续变量粗粒化，然后利用每个格子代替k值进行上述操作；
但是Logarithmic Binning中，横坐标取的是线性区间的中间值，而画图时为对数区间，所以这个上面可能需要有所商量，即有可能取对数区间最优；
特别值得注意的是，当我们把上述log-binning之后的区间画在\(\log\)图中的时候，其点的横坐标并非是等间隔均匀分布的，当格子的下标比较大的时候，横坐标才会逼近等间隔分布。用度分布的例子，证明如下： log-binning的格子为： \[ \begin{equation} k=1,k=2,3,k=4,5,6,7,...,k=2^{n-1},2^{n-1}+1,...,2^{n}-1 \end{equation} \] 平均度\(\left \langle k_{n} \right \rangle\)的分布为： \[ \begin{equation} \begin{aligned} \left \langle k_{n} \right \rangle&=\frac{S_{n-1}+1+S_{n}}{2}\\ &=\frac{2^{n-1}-1+1+2^{n}-1}{2}\\ &=3\times2^{n-2}-\frac{1}{2} \end{aligned} \end{equation} \] 当\(n\rightarrow\infty\)时，\(\left\langle k_{n} \right\rangle\)在\(\log\)坐标下的值为： \[ \begin{equation} \begin{aligned} \left\langle k_{n} \right\rangle_{log}&=\log\left( \left\langle k_{n} \right\rangle \right) \\ &=\log\left( 3\times2^{n-2}-\frac{1}{2} \right)\\ &\approx \log\left( 3\times2^{n-2} \right)\\ &\approx (n-2)\log\left( 2 \right)+ \log\left( 3\right)\\ &\approx n\log\left( 2 \right)+ \log\left( \frac{3}{4}\right)\\ &\approx n\log\left( 2 \right) \end{aligned} \end{equation} \]

即在\(\log\)坐标下，\(\left\langle k_{n} \right\rangle\)的间距为\(\log(2)\)。另外观察上述推导过程，最终的间距只与选择的比值\(q\)有关，即间距为\(\log(q)\)。下面对一般情况下的离散变量进行证明，假设初值\(b_{1}=a_{1}\)，比值\(q>1\)，则前\(n\)项的求和为： \[ \begin{equation} S_{n}=\frac{a_{1}(1-q^{n})}{1-q} \end{equation} \] 平均度\(\left \langle k_{n} \right \rangle\)的分布为： \[ \begin{equation} \begin{aligned} \left \langle k_{n} \right \rangle&=\frac{S_{n-1}+1+S_{n}}{2}\\ &=\frac{(1+q)a_{1}q^{n-1}}{2(q-1)}+\frac{a_{1}}{1-q}+\frac{1}{2}\\ &\approx \frac{(1+q)a_{1}q^{n-1}}{2(q-1)} \end{aligned} \end{equation} \]

当\(n\rightarrow\infty\)时，\(\left\langle k_{n} \right\rangle\)在\(\log\)坐标下的值为： \[ \begin{equation} \begin{aligned} \left\langle k_{n} \right\rangle_{log}&=\log\left( \left\langle k_{n} \right\rangle \right)\\ &=\log\left(\frac{(1+q)a_{1}q^{n-1}}{2(q-1)} \right)\\ &=\log(\frac{(q+1)a_{1}}{2(q-1)q})+n\log(q)\\ &\approx n\log(q) \end{aligned} \end{equation} \] 证毕。

附录

%计算log-binning
function [edges_exponent,N_hist_exponent]=LogBinning(N_hist,edges,first_bin,exponent_base)
%N_hist为数据，edges为横坐标（其length比N_hist多1），first_bin为横坐标第一个格子大小；
%exponent_base为比值；
%edges_exponent为横坐标，N_hist_exponent为纵坐标值；

first_bin=1;
exponent_base=2;

upper_limit=first_bin;    %代表总的格子数；
length_edges=length(edges);
count_exponent=1;
edges_exponent(1)=(edges(upper_limit+1)+edges(1))/2;
N_hist_exponent(1)=sum(N_hist(1:upper_limit))/upper_limit;
while (upper_limit+exponent_base^count_exponent+1)<=length_edges
    upper_limit1=upper_limit+exponent_base^count_exponent;
    count_exponent=count_exponent+1;
    edges_exponent(count_exponent)=(edges(upper_limit1+1)+edges(upper_limit+1))/2;
    N_hist_exponent(count_exponent)=sum(N_hist((upper_limit+1):upper_limit1))/(upper_limit1-upper_limit);
    upper_limit=upper_limit1;
end

upper_limit1=upper_limit+exponent_base^count_exponent;
count_exponent=count_exponent+1;
edges_exponent(count_exponent)=(2*edges(upper_limit+1)+(upper_limit1-upper_limit)*(edges(3)-edges(2)))/2;
N_hist_exponent(count_exponent)=sum(N_hist((upper_limit+1):end))/(upper_limit1-upper_limit);
upper_limit=upper_limit1;

edges_exponent=edges_exponent';
N_hist_exponent=N_hist_exponent';
end

其中作为在\(\log\)横坐标下的edges_exponent值的间隔为：

log10(edges_exponent(2:end))-log10(edges_exponent(1:(end-1)))

ans =

    0.6021    0.3979    0.3424    0.3203    0.3104    0.3056    0.3033

由于其比值选的也是\(2\)，所以最终的间距逐渐趋近为\(\log(2)=0.30103\)。图取线性横坐标中值

注意

如果横坐标为连续变量，则图中的前几个点（从左数）会往上翘起来，并且大部分点相对理论值会偏右（或上）。这是因为，此时由于对连续变量进行hist的时候，会取区间的整个数值的平均值，所以当这个区间内的前后有数值上的单调性（如果左高右低），此值会被区间内左侧的值拉高，高于区间中间位置的值的大小。解决的办法是把hist的区间数取大，进而减小区间的大小，减小拉高的高度。但是值得注意的是，这个区间也不能取的太小，太小的话可能会因为涨落的原因，使数据点低垂下去，或者翘起来。所以需要根据数据的具体情况进行调整。其实这也一定程度上显露了区间位置被选在线性区间中间这个操作的缺点。当我们取横坐标为\(\log\)之后的值的平均值的时候，可以克服这个缺点：

图取对数横坐标中值另外上面点在横坐标上的分布其实也不均匀。下面对一般情况下的离散变量（连续变量可以通过粗粒化转化为离散变量）进行证明，假设初值\(b_{1}=a_{1}\)，比值\(q>1\)，则前\(n\)项的求和为： \[ \begin{equation} S_{n}=\frac{a_{1}(1-q^{n})}{1-q} \end{equation} \]

当\(n\rightarrow\infty\)，且\(q>1\)时，\(\left\langle k_{n} \right\rangle\)在\(\log\)坐标下的值为： \[ \begin{equation} \begin{aligned} \left\langle k_{n} \right\rangle_{log}&=\frac{\log(S_{n-1})+\log(S_{n})}{2}\\ &=\frac{\log(S_{n-1}S_{n})}{2}\\ &=\frac{\log(\frac{a_{1}^{2}(1-q^{n-1})(1-q^{n})}{(1-q)^{2}})}{2}\\ &\approx \frac{\log(\frac{a_{1}^{2}q^{2n}}{q(1-q)^{2}})}{2}\\ &\approx \frac{2n\log(q)+\log(\frac{a_{1}^{2}}{q(1-q)^{2}})}{2}\\ &\approx n\log(q)+\frac{1}{2}\log(\frac{a_{1}^{2}}{q(1-q)^{2}})\\ &\approx n\log(q)\\ \end{aligned} \end{equation} \] 比较之前得到的在线性区间下的结果，当\(n\rightarrow\infty\)时，两者一致。另外当\(a_{1}\)较小，\(q\)较大时，点在横坐标上的分布会比较均匀。

代码为：

%计算log-binning，输出线性坐标下的横坐标和纵坐标；
function [edges_exponent,N_hist_exponent]=LogBinning(N_hist,edges,first_bin,exponent_base)
%N_hist为数据，edges为横坐标（其length比N_hist多1），first_bin为横坐标第一个格子大小；
%exponent_base为比值；
%edges_exponent为横坐标，N_hist_exponent为纵坐标值；

first_bin=1;
exponent_base=2;

upper_limit=first_bin;    %代表总的格子数；
length_edges=length(edges);
count_exponent=1;
edges_exponent(1)=10^((log10(edges(upper_limit+1))+log10(edges(1)))/2);
N_hist_exponent(1)=sum(N_hist(1:upper_limit))/upper_limit;
while (upper_limit+exponent_base^count_exponent+1)<=length_edges
    upper_limit1=upper_limit+exponent_base^count_exponent;
    count_exponent=count_exponent+1;
    edges_exponent(count_exponent)=10^((log10(edges(upper_limit1+1))+log10(edges(upper_limit+1)))/2);    %注意横坐标需要取对数坐标下的中间值，求完之后为了防止混淆，再退化为线性坐标；
    N_hist_exponent(count_exponent)=sum(N_hist((upper_limit+1):upper_limit1))/(upper_limit1-upper_limit);
    upper_limit=upper_limit1;
end

% 注意此时超出范围，需要额外处理；
upper_limit1=upper_limit+exponent_base^count_exponent;
count_exponent=count_exponent+1;
edges_exponent(count_exponent)=10^((log10(edges(upper_limit+1))+log10(edges(upper_limit+1)+(upper_limit1-upper_limit)*(edges(3)-edges(2))))/2);
N_hist_exponent(count_exponent)=sum(N_hist((upper_limit+1):end))/(upper_limit1-upper_limit);
upper_limit=upper_limit1;

edges_exponent=edges_exponent';
N_hist_exponent=N_hist_exponent';
end

参考文献：barabasi，network science，chapter4.