使用NumPy进行 y = x ** a 的回归拟合

Question

import numpy as np
import matplotlib.pyplot as plt
x_ = np.array([1.0, 2.0, 4.0, 8.0, 16.0, 33.0])
y_ = np.array([0.4, 0.55, 0.62, 0.72, 0.76, 0.8])

我希望拟合一个回归模型，形式为y = x ** a，其中估计出a的值，目的是对(更)大的未知x值进行外推。注意: 我需要强制我的模型经过(0, 0)点。

不确定这是否是最好的方法，但我尝试将其转换为找到最适合log(y) / log(x)的常数。我已经:

x_ = np.array([1.0, 2.0, 4.0, 8.0, 16.0, 33.0])
logx_ = np.log(x_)
y_ = np.array([0.4, 0.55, 0.62, 0.72, 0.76, 0.8])
logy_ = np.log(y_)
y = logy_ / logx_
x = x_
print "y: ", y
z = np.polyfit(x, y, 0)
print "param estimates: ", z

p = np.poly1d(z)
p30 = np.poly1d(np.polyfit(x, y, 30))

xp = np.linspace(0, 35, 6)
_ = plt.plot(x, np.exp(y*np.log(x)), '.', xp, np.exp(p(xp)*np.log(x)), '-')
plt.ylim(0,1.0)
plt.show()

但似乎并没有起作用。 有人能解释一下我做错了什么，并提供一个答案吗？

- Alexandre Holden Daly

---

最好将线性拟合到ln(y)而不是ln(x)，因为根据您的模型ln(y)=a*ln(x)。然后你需要斜率。但我不确定你能否强制polyfit跳过一个常数；你可能需要使用scipy.optimize.curve_fit（但这样你可以直接拟合模型）。 - Andras Deak -- Слава Україні

Answer 1

使用:

log(y+1) = a * log(x+1)

这将帮助您获得指数函数并通过点 (0, 0)。

f = lambda x, a: (x + 1) ** a - 1
x = np.random.rand(200, 1) * 100
y = f(x, 1.5) + np.random.rand(*x.shape) * 200

plt.plot(x, y, '.')

log_x_plus_1 = np.log(x + 1)
log_y_plus_1 = np.log(y + 1)
a = np.linalg.pinv(log_x_plus_1.T.dot(log_x_plus_1)).dot(log_x_plus_1.T).dot(log_y_plus_1)[0][0]

plt.plot(x, y, '.')
plt.plot(x, f(x, a), 'r.')

- piRSquared

---

你的意思是 log(y+1) = a * log(x+1) 吗？ - Alicia Garcia-Raboso

是的！我会的。我有机会时会更新帖子。这个想法在那里。 - piRSquared

但在这种情况下，你不是在拟合 y = x ** a 而是在拟合 y+1 = (x+1) ** a... - Alicia Garcia-Raboso

Answer 2

通过对 logy_ = a * logx_ 进行最小二乘法拟合，可以得到一个纯的 y = x ** a 模型：

a = np.linalg.lstsq(logx_[:, np.newaxis], logy_)[0][0]

- Alicia Garcia-Raboso