我有一组数据，我想比较哪条线最能描述它(不同阶的多项式，指数或对数)。
我使用Python和Numpy和多项式拟合有一个函数polyfit()。但我没有发现这样的函数的指数和对数拟合。
有吗?或者如何解决这个问题?
当前回答
下面是使用scikit learn中的工具对简单数据进行线性化的选项。
鉴于
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import FunctionTransformer
np.random.seed(123)
# General Functions
def func_exp(x, a, b, c):
"""Return values from a general exponential function."""
return a * np.exp(b * x) + c
def func_log(x, a, b, c):
"""Return values from a general log function."""
return a * np.log(b * x) + c
# Helper
def generate_data(func, *args, jitter=0):
"""Return a tuple of arrays with random data along a general function."""
xs = np.linspace(1, 5, 50)
ys = func(xs, *args)
noise = jitter * np.random.normal(size=len(xs)) + jitter
xs = xs.reshape(-1, 1) # xs[:, np.newaxis]
ys = (ys + noise).reshape(-1, 1)
return xs, ys
transformer = FunctionTransformer(np.log, validate=True)
Code
拟合指数数据
# Data
x_samp, y_samp = generate_data(func_exp, 2.5, 1.2, 0.7, jitter=3)
y_trans = transformer.fit_transform(y_samp) # 1
# Regression
regressor = LinearRegression()
results = regressor.fit(x_samp, y_trans) # 2
model = results.predict
y_fit = model(x_samp)
# Visualization
plt.scatter(x_samp, y_samp)
plt.plot(x_samp, np.exp(y_fit), "k--", label="Fit") # 3
plt.title("Exponential Fit")
拟合日志数据
# Data
x_samp, y_samp = generate_data(func_log, 2.5, 1.2, 0.7, jitter=0.15)
x_trans = transformer.fit_transform(x_samp) # 1
# Regression
regressor = LinearRegression()
results = regressor.fit(x_trans, y_samp) # 2
model = results.predict
y_fit = model(x_trans)
# Visualization
plt.scatter(x_samp, y_samp)
plt.plot(x_samp, y_fit, "k--", label="Fit") # 3
plt.title("Logarithmic Fit")
细节
一般步骤
对数据值(x、y或两者)应用日志操作 将数据回归到线性化模型 通过“反转”任何日志操作(使用np.exp())绘制并适合原始数据
假设我们的数据遵循指数趋势，一般的方程+可能是:
我们可以通过取对数来线性化后一个方程(例如y =截距+斜率* x):
给定线性化方程++和回归参数，我们可以计算出:
A经截距(ln(A)) B经斜率(B)
线性化技术概述
Relationship | Example | General Eqn. | Altered Var. | Linearized Eqn.
-------------|------------|----------------------|----------------|------------------------------------------
Linear | x | y = B * x + C | - | y = C + B * x
Logarithmic | log(x) | y = A * log(B*x) + C | log(x) | y = C + A * (log(B) + log(x))
Exponential | 2**x, e**x | y = A * exp(B*x) + C | log(y) | log(y-C) = log(A) + B * x
Power | x**2 | y = B * x**N + C | log(x), log(y) | log(y-C) = log(B) + N * log(x)
+注意:线性化指数函数在噪声较小且C=0时效果最好。请谨慎使用。
注:修改x数据有助于线性化指数数据，而修改y数据有助于线性化日志数据。
其他回答
我们在解决这两个问题的同时演示了lmfit的特性。
鉴于
import lmfit
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
np.random.seed(123)
# General Functions
def func_log(x, a, b, c):
"""Return values from a general log function."""
return a * np.log(b * x) + c
# Data
x_samp = np.linspace(1, 5, 50)
_noise = np.random.normal(size=len(x_samp), scale=0.06)
y_samp = 2.5 * np.exp(1.2 * x_samp) + 0.7 + _noise
y_samp2 = 2.5 * np.log(1.2 * x_samp) + 0.7 + _noise
Code
方法1 - lmfit模型
拟合指数数据
regressor = lmfit.models.ExponentialModel() # 1
initial_guess = dict(amplitude=1, decay=-1) # 2
results = regressor.fit(y_samp, x=x_samp, **initial_guess)
y_fit = results.best_fit
plt.plot(x_samp, y_samp, "o", label="Data")
plt.plot(x_samp, y_fit, "k--", label="Fit")
plt.legend()
方法2 -自定义模型
拟合日志数据
regressor = lmfit.Model(func_log) # 1
initial_guess = dict(a=1, b=.1, c=.1) # 2
results = regressor.fit(y_samp2, x=x_samp, **initial_guess)
y_fit = results.best_fit
plt.plot(x_samp, y_samp2, "o", label="Data")
plt.plot(x_samp, y_fit, "k--", label="Fit")
plt.legend()
细节
选择一个回归类 提供命名的、尊重函数域的初始猜测
您可以从回归器对象确定推断出的参数。例子:
regressor.param_names
# ['decay', 'amplitude']
要进行预测，请使用ModelResult.eval()方法。
model = results.eval
y_pred = model(x=np.array([1.5]))
注意:ExponentialModel()遵循一个衰减函数，它接受两个参数，其中一个是负的。
请参见ExponentialGaussianModel()，它接受更多参数。
通过> pip Install lmfit安装库。
下面是使用scikit learn中的工具对简单数据进行线性化的选项。
鉴于
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import FunctionTransformer
np.random.seed(123)
# General Functions
def func_exp(x, a, b, c):
"""Return values from a general exponential function."""
return a * np.exp(b * x) + c
def func_log(x, a, b, c):
"""Return values from a general log function."""
return a * np.log(b * x) + c
# Helper
def generate_data(func, *args, jitter=0):
"""Return a tuple of arrays with random data along a general function."""
xs = np.linspace(1, 5, 50)
ys = func(xs, *args)
noise = jitter * np.random.normal(size=len(xs)) + jitter
xs = xs.reshape(-1, 1) # xs[:, np.newaxis]
ys = (ys + noise).reshape(-1, 1)
return xs, ys
transformer = FunctionTransformer(np.log, validate=True)
Code
拟合指数数据
# Data
x_samp, y_samp = generate_data(func_exp, 2.5, 1.2, 0.7, jitter=3)
y_trans = transformer.fit_transform(y_samp) # 1
# Regression
regressor = LinearRegression()
results = regressor.fit(x_samp, y_trans) # 2
model = results.predict
y_fit = model(x_samp)
# Visualization
plt.scatter(x_samp, y_samp)
plt.plot(x_samp, np.exp(y_fit), "k--", label="Fit") # 3
plt.title("Exponential Fit")
拟合日志数据
# Data
x_samp, y_samp = generate_data(func_log, 2.5, 1.2, 0.7, jitter=0.15)
x_trans = transformer.fit_transform(x_samp) # 1
# Regression
regressor = LinearRegression()
results = regressor.fit(x_trans, y_samp) # 2
model = results.predict
y_fit = model(x_trans)
# Visualization
plt.scatter(x_samp, y_samp)
plt.plot(x_samp, y_fit, "k--", label="Fit") # 3
plt.title("Logarithmic Fit")
细节
一般步骤
对数据值(x、y或两者)应用日志操作 将数据回归到线性化模型 通过“反转”任何日志操作(使用np.exp())绘制并适合原始数据
假设我们的数据遵循指数趋势，一般的方程+可能是:
我们可以通过取对数来线性化后一个方程(例如y =截距+斜率* x):
给定线性化方程++和回归参数，我们可以计算出:
A经截距(ln(A)) B经斜率(B)
线性化技术概述
Relationship | Example | General Eqn. | Altered Var. | Linearized Eqn.
-------------|------------|----------------------|----------------|------------------------------------------
Linear | x | y = B * x + C | - | y = C + B * x
Logarithmic | log(x) | y = A * log(B*x) + C | log(x) | y = C + A * (log(B) + log(x))
Exponential | 2**x, e**x | y = A * exp(B*x) + C | log(y) | log(y-C) = log(A) + B * x
Power | x**2 | y = B * x**N + C | log(x), log(y) | log(y-C) = log(B) + N * log(x)
+注意:线性化指数函数在噪声较小且C=0时效果最好。请谨慎使用。
注:修改x数据有助于线性化指数数据，而修改y数据有助于线性化日志数据。
Wolfram有一个封闭形式的解拟合指数。他们也有类似的解决方案拟合对数和幂律。
我发现这比scipy的curve_fit工作得更好。特别是当你没有“接近零”的数据时。这里有一个例子:
import numpy as np
import matplotlib.pyplot as plt
# Fit the function y = A * exp(B * x) to the data
# returns (A, B)
# From: https://mathworld.wolfram.com/LeastSquaresFittingExponential.html
def fit_exp(xs, ys):
S_x2_y = 0.0
S_y_lny = 0.0
S_x_y = 0.0
S_x_y_lny = 0.0
S_y = 0.0
for (x,y) in zip(xs, ys):
S_x2_y += x * x * y
S_y_lny += y * np.log(y)
S_x_y += x * y
S_x_y_lny += x * y * np.log(y)
S_y += y
#end
a = (S_x2_y * S_y_lny - S_x_y * S_x_y_lny) / (S_y * S_x2_y - S_x_y * S_x_y)
b = (S_y * S_x_y_lny - S_x_y * S_y_lny) / (S_y * S_x2_y - S_x_y * S_x_y)
return (np.exp(a), b)
xs = [33, 34, 35, 36, 37, 38, 39, 40, 41, 42]
ys = [3187, 3545, 4045, 4447, 4872, 5660, 5983, 6254, 6681, 7206]
(A, B) = fit_exp(xs, ys)
plt.figure()
plt.plot(xs, ys, 'o-', label='Raw Data')
plt.plot(xs, [A * np.exp(B *x) for x in xs], 'o-', label='Fit')
plt.title('Exponential Fit Test')
plt.xlabel('X')
plt.ylabel('Y')
plt.legend(loc='best')
plt.tight_layout()
plt.show()
还可以使用scipy.optimize中的curve_fit将一组数据适合于任何函数。例如，如果你想拟合一个指数函数(来自文档):
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
def func(x, a, b, c):
return a * np.exp(-b * x) + c
x = np.linspace(0,4,50)
y = func(x, 2.5, 1.3, 0.5)
yn = y + 0.2*np.random.normal(size=len(x))
popt, pcov = curve_fit(func, x, yn)
然后如果你想画出来，你可以这样做:
plt.figure()
plt.plot(x, yn, 'ko', label="Original Noised Data")
plt.plot(x, func(x, *popt), 'r-', label="Fitted Curve")
plt.legend()
plt.show()
(注意:当你绘图时，popt前面的*会将项展开为func所期望的a, b和c。)
我想你可以用:
np.log --> natural log
np.log10 --> base 10
np.log2 --> base 2
稍微修改一下IanVS的答案:
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
def func(x, a, b, c):
#return a * np.exp(-b * x) + c
return a * np.log(b * x) + c
x = np.linspace(1,5,50) # changed boundary conditions to avoid division by 0
y = func(x, 2.5, 1.3, 0.5)
yn = y + 0.2*np.random.normal(size=len(x))
popt, pcov = curve_fit(func, x, yn)
plt.figure()
plt.plot(x, yn, 'ko', label="Original Noised Data")
plt.plot(x, func(x, *popt), 'r-', label="Fitted Curve")
plt.legend()
plt.show()
结果如下图所示:
推荐文章
- 美丽的汤和提取div及其内容的ID
- 在Python中重置生成器对象
- 用Python构建最小的插件架构
- model.eval()在pytorch中做什么?
- Tensorflow 2.0:模块“Tensorflow”没有属性“Session”
- 从环境文件中读入环境变量
- 在OSX 10.11中安装Scrapy时，“OSError: [Errno 1]操作不允许”(El Capitan)(系统完整性保护)
- 如何删除熊猫数据帧的最后一行数据
- 我如何在熊猫中找到数字列?
- 检查pandas数据框架索引中是否存在值
- 计算熊猫数量的最有效方法是什么?
- 如何在python中验证日期字符串格式?
- 用csv模块从csv文件中读取特定的列?
- 使用PyCrypto AES-256加密和解密
- “\d”在正则表达式中是数字吗?
