布朗运动可见度图的幂律拟合优化方案问询

尝试复现论文《From time series to complex networks: The visibility graph》的分析流程：生成布朗运动B(t)数据，通过ts2vg构建自然可见度图，计算度分布后分别用powerlaw库和polyfit函数做幂律拟合。但度分布存在概率为0的取值，调整后拟合效果仍不理想，寻求优化建议及可达到论文拟合效果的替代方法。

数据生成代码 #

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
def brownian_motion(t, a=1):
    dt = t[1] - t[0]
    dW = np.random.normal(scale=np.sqrt(dt), size=t.shape)
    W = np.cumsum(dW)
    return a**(1/2) * W

t = np.linspace(0, 100, 100000)

np.random.seed(42)

W = brownian_motion(t)

plt.plot(t, W)
plt.xlabel('t')
plt.ylabel('W(t)')
plt.title('Brownian Motion')
plt.show()

df = pd.DataFrame(W, index=t, columns=['W(t)'])
df.index.name = 't'
df.to_csv('brownian_motion_data.csv')

拟合代码 #

import numpy as np
import matplotlib.pyplot as plt
from ts2vg import NaturalVG
import networkx as nx
import pandas as pd
import powerlaw

data = pd.read_csv('brownian_motion_data.csv')
dat = data['W(t)'][:500]
nodes = range(len(dat))
vg = NaturalVG()
vg.build(dat)
edges = vg.edges
G = nx.Graph()
G.add_nodes_from(nodes)
G.add_edges_from(edges)

n = nx.number_of_nodes(G)
degree_freq = nx.degree_histogram(G)
degree_dist = []
for i in degree_freq:
   x = i / n
   degree_dist.append(x)
degrees = range(len(degree_freq))
fit = powerlaw.Fit(np.array(degree_dist) + 10 ** (-6))
alpha = round(fit.power_law.alpha, 3)

plt.figure(figsize=(8, 6))
plt.scatter(degrees, degree_dist, marker='o', color='b', label='Data')
fit.power_law.plot_pdf(color='r', linestyle='--', label=f'Fit (alpha={alpha})')
plt.xlabel('Degree')
plt.ylabel('Frequency')
plt.legend()
plt.show()

degrees_array = np.array(degrees)
degree_dist_array = np.array(degree_dist)
coefficients = np.polyfit(degrees_array, degree_dist_array, 2)
poly_function = np.poly1d(coefficients)
degrees_fit = np.linspace(min(degrees_array), max(degrees_array), 100)
degree_dist_fit = poly_function(degrees_fit)

plt.loglog(degrees_array, degree_dist_array, label='Original Data')
plt.loglog(degrees_fit, degree_dist_fit, label=f'Polynomial Fit', color='red')
plt.xlabel('Degrees')
plt.ylabel('Degree Distribution')
plt.title('Polynomial Fitting Example for Degree Distribution')
plt.legend()
plt.show()

一、度分布预处理优化 #

• 过滤零概率度值：直接剔除度分布中概率为0的点，只保留有实际出现频率的(度值, 概率)对，避免无效数据干扰拟合。示例代码：

  # 过滤零概率点
  valid_degrees = []
  valid_probs = []
  for d, p in zip(degrees, degree_dist):
      if p > 0:
          valid_degrees.append(d)
          valid_probs.append(p)
  valid_degrees = np.array(valid_degrees)
  valid_probs = np.array(valid_probs)
  

• 使用原始计数而非概率：powerlaw库更适合输入度的原始出现次数（即degree_freq中大于0的部分），而非归一化后的概率，避免人为添加小扰动带来的偏差。

二、powerlaw库拟合参数调整 #

• 指定拟合范围：通过xmin参数手动指定幂律拟合的起始度值，或让库自动选择最优xmin，聚焦在符合幂律特征的区间：

  # valid_counts为过滤后的度计数数组
  fit = powerlaw.Fit(valid_counts, xmin='estimate')
  

• 对比模型拟合效果：用fit.distribution_compare()方法对比幂律模型与指数、对数正态等模型的拟合优度，确认数据是否真的符合幂律，避免强行拟合。

三、替代拟合方法 #

• 手动实现最大似然估计(MLE)：针对幂律分布$P(k) \propto k^{-\alpha}$，直接通过MLE计算指数：

  # k_values为过滤后的有效度值数组
  alpha = 1 + len(k_values) / np.sum(np.log(k_values / min(k_values)))
  

• 双对数坐标下线性拟合：幂律分布在双对数图中呈直线，对过滤后的双对数数据做线性拟合（而非二次多项式）：

  log_k = np.log10(valid_degrees)
  log_p = np.log10(valid_probs)
  coeffs = np.polyfit(log_k, log_p, 1)
  alpha = -coeffs[0]  # 幂律指数为斜率绝对值
  fit_line = np.poly1d(coeffs)(log_k)
  
  plt.loglog(valid_degrees, valid_probs, 'bo', label='Data')
  plt.loglog(valid_degrees, 10**fit_line, 'r--', label=f'Fit (alpha={round(alpha,3)})')
  plt.legend()
  plt.show()
  

• NetworkX内置拟合工具：用nx.degree_power_law_fit()直接对图的度分布做幂律拟合，简化流程：

  alpha, xmin = nx.degree_power_law_fit(G)
  

四、数据与模型构建优化 #

• 增加序列长度：当前仅使用前500个数据点，论文可能采用更长序列，增加数据量能让度分布更稳定，接近理论幂律特征。
• 验证可见度图正确性：检查ts2vg生成的边是否符合自然可见度规则：任意节点$i,j$可见当且仅当中间所有节点$k$满足$W(k) < W(i) + (W(j)-W(i))*(k-i)/(j-i)$，可手动抽查小范围节点连接情况。

内容的提问来源于stack exchange，提问作者Dinh Duc Vu