Python2.7中日志日志刻度的最佳拟合

Question

这是一个网络IP频率等级图，以日志为单位。在完成这一部分之后，我将尝试使用python2.7在日志刻度上绘制最佳拟合线。我必须使用matplotlib的“symlog”轴比例，否则某些值将无法正确显示，某些值将被隐藏。在

我绘制的数据的X值是url，Y值是url的相应频率。在

我的数据如下：

'http://www.bing.com/search?q=d2l&src=IE-TopResult&FORM=IETR02&conversationid=  123 0.00052210688591'
`http://library.uc.ca/  118 4.57782298326e-05`
`http://www.bing.com/search?q=d2l+uofc&src=IE-TopResult&FORM=IETR02&conversationid= 114 4.30271029472e-06`
`http://www.nature.com/scitable/topicpage/genetics-and-statistical-analysis-34592   109 1.9483268261e-06`

数据在第一列中包含URL，在第二列中包含相应的频率（同一URL出现的次数），最后在第三列中包含传输的字节。首先，我只使用第一列和第二列进行分析。总共有2465x值或唯一的url。在

下面是我的代码

^{pr2}$

你会看到大量的库被导入，因为我一直在使用它们，但我的实验都没有产生预期的结果。所以上面的代码正确地生成了秩图。这是红线，但曲线中的蓝线，被认为是最佳拟合线，在视觉上是不正确的，可以看出。这是生成的图形。在

这是我期待的图表。第二张图中的虚线是我不正确地绘制的。在

有什么办法解决这个问题吗？在

Answer 1

在日志刻度上沿直线下降的数据遵循y = c*x^(m)形式的幂关系。通过取两边的对数，可以得到要拟合的线性方程：

log(y) = m*log(x) + c

调用np.polyfit(log(x), log(y), 1)提供m和{}的值。然后可以使用这些值计算log_y_fit的拟合值，如下所示：

^{pr2}$

要根据原始数据绘制的拟合值为：

y_fit = exp(log_y_fit) = exp(m*log(x) + c)

所以，你面临的两个问题是：

1. 使用原始x坐标而不是对数（x）坐标计算拟合值

2. 绘制拟合y值的对数而不将其转换回原始比例

在下面的代码中，我将plt.plot(z, np.poly1d(np.polyfit(logA, logB, 1))(z))替换为：

m, c = np.polyfit(logA, logB, 1) # fit log(y) = m*log(x) + c
y_fit = np.exp(m*logA + c) # calculate the fitted values of y 
plt.plot(z, y_fit, ':')

这可以放在一行：plt.plot(z, np.exp(np.poly1d(np.polyfit(logA, logB, 1))(logA)))，但我发现这使得调试更加困难。在

以下代码中有一些不同之处：

• 当您从z计算logA以过滤掉任何值<；1时，您正在使用列表理解，但是z是一个线性范围，只有第一个值是<；1。从1开始创建z似乎更容易，我就是这样编码的。

• 我不知道你为什么在你的列表理解中使用x*log(x)这个词。这在我看来是个错误，所以我没有把它包括在答案中。

此代码应能正确工作：

fig=plt.figure()
ax = fig.add_subplot(111)

z=np.arange(1, len(x)+1) #start at 1, to avoid error from log(0)

logA = np.log(z) #no need for list comprehension since all z values >= 1
logB = np.log(y)

m, c = np.polyfit(logA, logB, 1) # fit log(y) = m*log(x) + c
y_fit = np.exp(m*logA + c) # calculate the fitted values of y 

plt.plot(z, y, color = 'r')
plt.plot(z, y_fit, ':')

ax.set_yscale('symlog')
ax.set_xscale('symlog')
#slope, intercept = np.polyfit(logA, logB, 1)
plt.xlabel("Pre_referer")
plt.ylabel("Popularity")
ax.set_title('Pre Referral URL Popularity distribution')
plt.show()

当我在模拟数据上运行它时，我得到以下图表：

注意事项：

• 行的左右两端的“扭结”是使用“symlog”的结果，它将极小的值线性化，如What is the difference between 'log' and 'symlog'?的答案所述。如果这些数据绘制在“对数-对数”轴上，则拟合的数据将是一条直线。

• 您可能还想阅读这个答案：https://stackoverflow.com/a/3433503/7517724，它解释了如何使用加权来实现对日志转换数据的“更好”拟合。

Answer 2

我想出了解决这个问题的另一个办法。分享这个是因为它可能有帮助。在

fig=plt.figure()
ax = fig.add_subplot(111)
z=np.arange(len(x)) + 1
print z
print y
rank = [np.log10(i) for i in z]
freq = [np.log10(i) for i in y]
m, b, r_value, p_value, std_err = stats.linregress(rank, freq)
print "slope: ", m
print "r-squared: ", r_value**2
print "intercept:", b
plt.plot(rank, freq, 'o',color = 'r')
abline_values = [m * i + b for i in rank]
plt.plot(rank, abline_values)

这基本上也达到了目标。它使用stats模块。在