python 2.7中对数刻度的最佳拟合线

Posted 2023-03-12

【中文标题】python 2.7中对数刻度的最佳拟合线

【英文标题】：Best Fit Line on Log Log Scales in python 2.7

【发布时间】：2017-10-05 20:07:01

【问题描述】：

这是对数刻度的网络 IP 频率等级图。完成这部分后，我尝试使用 Python 2.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 出现的次数），最后在第三列中包含传输的字节数。首先，我只使用第 1 列和第 2 列进行此分析。共有 2,465 个 x 值或唯一 URL。

以下是我的代码

import os
import matplotlib.pyplot as plt
import numpy as np
import math
from numpy import *
import scipy
from scipy.interpolate import *
from scipy.stats import linregress
from scipy.optimize import curve_fit

file = open(filename1, 'r')
lines = file.readlines()

result = 
x=[]
y=[]
for line in lines:
  course,count,size = line.lstrip().rstrip('\n').split('\t')
  if course not in result:
      result[course] = int(count)
  else:
      result[course] += int(count)
file.close()

frequency = sorted(result.items(), key = lambda i: i[1], reverse= True)
x=[]
y=[]
i=0
for element in frequency:
  x.append(element[0])
  y.append(element[1])


z=[]
fig=plt.figure()
ax = fig.add_subplot(111)
z=np.arange(len(x))
print z
logA = [x*np.log(x) if x>=1 else 1 for x in z]
logB = np.log(y)
plt.plot(z, y, color = 'r')
plt.plot(z, np.poly1d(np.polyfit(logA, logB, 1))(z))
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()

您会看到导入了很多库，因为我一直在使用其中的很多库，但我的实验都没有产生预期的结果。所以上面的代码正确地生成了排名图。可以看出，这是红线，但曲线中应该是最佳拟合线的蓝线在视觉上是不正确的。这是生成的图表。

这是我期待的图表。第二张图中的虚线是我以某种方式错误地绘制的。

关于如何解决此问题的任何想法？

【参考方案1】：

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

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

调用np.polyfit(log(x), log(y), 1) 提供m 和c 的值。然后，您可以使用这些值来计算 log_y_fit 的拟合值：

log_y_fit = m*log(x) + c

您想要根据原始数据绘制的拟合值是：

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

所以，您遇到的两个问题是：

您正在使用原始 x 坐标计算拟合值，而不是 log(x) 坐标

您正在绘制拟合 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 以过滤掉任何 z 是一个线性范围，只有第一个值是 z，这就是我的编码方式。

我不知道你为什么在logA 的列表理解中有术语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'? 的答案中所述。如果此数据绘制在“log-log”轴上，则拟合数据将是一条直线。

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

【参考方案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)

这也基本上实现了目标。它使用统计模块。