pandas Python将多项式拟合到自然对数，然后将转换反转回非对数空间？

6rvt4ljy 于 11个月前发布在 Python

关注(0)|答案(2)|浏览(71)

我试图创建一个平滑函数来表示我正在处理的一些数据。问题是数据非常嘈杂，简单地使用最小二乘法来创建最佳拟合的多项式线并不那么好。

polynomial = np.polyfit(df.x,df.y)

字符串
做一些挖掘，显然对于这类事情，有些人首先将他们的输入数据转换为自然对数，这会产生一个更好的曲线。

df['y']=np.log(df['y'])
lop_polynomial = np.polyfit(df.x,df.y)

型
当我把这个多项式转换成一个函数并画出一堆点（只取输出的e^x）时，它看起来就像我想要的那样;

log_polynomial =np.poly1d(put_log_polynomial)

x_values = np.arange(0, 100,.01)
plt.plot(
    x_values,
    np.exp(log_polynomial (x_values ))

型
然而，这需要我将输入数据转换为日志数据，然后找到一个最佳拟合的多项式，然后将该多项式转换为函数，然后将该函数转换为一组离散点，然后在这些点周围创建另一个最佳拟合的多项式。因为这是我要运行很多次的事情，这是很多额外的步骤。
有没有一种方法可以将我使用对数输入得到的多项式转换回非对数空间，而无需手动对所有输入进行指数运算以获得离散点，然后尝试为非对数空间中的这些点创建一条最佳拟合线？

**edit * 这里是我的用例的更多背景。我有一组离散点，它们形成一条曲线。这条曲线的二阶导数包含了我想要干净地提取的重要信息，但是如果我不加修改地取二阶导数，它最终会非常嘈杂，将曲线中的小扭结过度解释为二阶导数中的巨大粗糙尖峰。

我尝试了各种方法来平滑曲线，但我尝试的所有其他方法要么无法充分解决噪声问题，要么破坏了我试图隔离的固有信号。以下是我迄今为止尝试的方法：
非对数空间中最佳拟合的最小二乘直线

polynomial = np.polyfit(df.x,df.y)

型
一维高斯滤波器

scipy.ndimage.gaussian_filter1d(df.y, 3)

型
三次样条

t,c,k=scipy.interpolate.splrep(df.x,df.y,k=3)

型
B-spline

spline = scipy.interpolate.BSpline(t, c, k, extrapolate=False)

型
我在找到最佳拟合多项式之前尝试对其进行对数变换的原因是，我发现一些研究论文使用这种插值方法来处理我正在研究的确切用例，当我在我的数据上测试它时，它似乎工作得很好。

edit 2**正如@NickODell敏锐地指出的那样，有很多潜在的多项式不能转换回非对数空间，仍然以python多项式格式表示x的整数幂的系数列表。因此，据我所知，这使得将对数空间多项式直接有效地转换回非对数空间中的多项式是不可能的。

我真正想要的是输入数据的二阶导数，它没有噪声，保留了所有的信号。我希望我可以通过将对数空间多项式直接转换回非对数空间并直接导出它来做到这一点，但显然这是不可能的。然而，它仍然可以表示为一个函数，我可以使用np.gradient导出，这将有效地允许同样的事情。但我不太相信这实际上会显着提高程序的效率。我将尝试弄乱一些分析，看看我是否能得到一个处理多大的交易。

edit 3**我要感谢@jlandercy和@Nick奥德尔清晰透彻的回答，我希望我能选择他们两个作为答案。

pandas

来源：https://stackoverflow.com/questions/77708524/python-fitting-polynomial-to-natural-log-then-reversing-the-transformation-back

2条答案

按热度按时间

s8vozzvw1#

如果我正确理解了你的要求，你想要：

在y轴上应用对数变换;
拟合噪声较大的多项式数据;
预测新x点的值;
对拟合曲线求二阶导数;
使必须多次执行的过程自动化。

Sklearn是一个通过Pipeline s执行此类操作的完美软件包。
首先，我们创建一些多项式噪声数据：

import numpy as np

from sklearn.preprocessing import PolynomialFeatures
from sklearn.compose import TransformedTargetRegressor
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline

np.random.seed(12345)
x = np.linspace(-1, 1, 50).reshape(-1, 1)
P = np.poly1d([6, 5, 4, 3, 2, 10])
y = P(x[:, 0])
n = 0.5 * np.random.randn(y.size)
yn = y + n

字符串

简单回归

作为基线，我们创建一个管道，从您的数据中线性回归多项式：

pipeline_lin = Pipeline([
    ("transformer", PolynomialFeatures(5)),
    ("model", LinearRegression(fit_intercept=False))
])

型
并适合数据集：

pipeline_lin.fit(x, y)

pipeline_lin["model"].coef_
# array([10.0226118 ,  1.51830909,  2.23638955,  2.91665499,  5.99904944, 7.88614861])

型

日志转换

如果你想在目标上应用对数变换（y轴，实际上它需要正值），我们可以改变管道：

pipeline_log = Pipeline([
    ("transformer", PolynomialFeatures(5)),
    ("model",
         TransformedTargetRegressor(
             regressor=LinearRegression(fit_intercept=False),
             func=np.log,
             inverse_func=np.exp
         )
    )
])
pipeline_log.fit(x, yn)

pipeline_log["model"].regressor_.coef_
# array([ 2.2901081 ,  0.15049063,  0.46685674,  0.32678866, -0.12522207, 0.34130133])

型
实际上，系数与原始多项式不直接相关，因为我们应用了对数变换。
现在我们可以预测其他特征的新目标值，而无需关心前向和后向日志转换：

xlin = np.linspace(-1, 1, 200).reshape(-1, 1)
yh_lin = pipeline_lin.predict(xlin)
yh_log = pipeline_log.predict(xlin)

fig, axe = plt.subplots()
axe.scatter(x, yn, marker=".", label="Data")
axe.plot(xlin, yh_lin, label="Linear")
axe.plot(xlin, yh_log, label="Log Linear")
axe.legend()
axe.grid()

型

的数据

自动化

这个过程可以重复多次，只需链接fit和predict：

# Iterate datasets:
for i in range(10):
    x = ...   # Select X data
    yn = ...  # Select Y data
    # Fit and predict:
    yhat = pipeline_log.fit(x, yn).predict(xlin)

型

可选

如果你不知道多项式的阶数，也不要求拟合是多项式回归（例如，你只想预测合理的值），那么高斯过程可以是一个很好的选择：

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF

pipeline_gpr = Pipeline([
    ("model",
         TransformedTargetRegressor(
             regressor=GaussianProcessRegressor(
                 kernel=1.0*RBF(),
                 alpha=0.01
             ),
             func=np.log,
             inverse_func=np.exp
         )
    )
])

型

的

求二阶导数

编号
估计二阶导数的一个潜在候选是从scipy得到savgol_filter：

from scipy import signal

dx = xlin[1,0] - xlin[0,0]
yh_lin_d2 = signal.savgol_filter(yh_lin, 15, 3, deriv=2, delta=dx)
yh_log_d2 = signal.savgol_filter(yh_log, 15, 3, deriv=2, delta=dx)
yh_gpr_d2 = signal.savgol_filter(yh_gpr, 15, 3, deriv=2, delta=dx)

型

注意，取拟合的导数（或者更糟糕的是二阶导数）不是一个简单的操作。正如预期的那样，5阶多项式（线性）的导数是3阶多项式。

根据回归系数

若拟合模型为多项式，则可利用回归多项式系数推导出任意阶导数。

Phat = np.poly1d(pipeline["model"].coef_[::-1])
Phat_xx = Phat.deriv(m=2)

型
可以用来插入新点：

yh_d2 = Phat_xx(np.squeeze(xlin))

型
这个解决方案完全符合yh_lin_d2。

提案

下面是一个完整的建议，基于一组遵循多项式趋势的点来求二阶导数：

def fit_diff(x, y, poly_order=5, resolution=200):
    pipeline = Pipeline([
        ("transformer", PolynomialFeatures(poly_order)),
        ("model", LinearRegression(fit_intercept=False))
    ])
    xlin = np.linspace(x.min(), x.max(), resolution).reshape(-1, 1)
    yhat = pipeline.fit(x, y).predict(xlin)
    dx = xlin[1, 0] - xlin[0, 0]
    dy2dx2 = signal.savgol_filter(yhat, 15, 3, deriv=2, delta=dx)
    return dy2dx2

型

赞(0）回复(0）举报 11个月前

xwbd5t1u2#

除了@jlandercy建议的数值解之外，还有一个解析解。不幸的是，正如我在评论中讨论的那样，它不是一个很好的多项式。因此，我们需要引入一个新的库SymPy，它能够表示，简化和区分任意表达式。
我将从创建他使用的相同数据集开始，以使答案更容易比较。

import sympy as sm
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(12345)
x_data = np.linspace(-1, 1, 50)
P = np.poly1d([6, 5, 4, 3, 2, 10])
y_data = P(x_data) + 0.5 * np.random.randn(x_data.size)

字符串
然后，将多项式拟合到记录的数据版本：

poly_coeffs = np.polyfit(x_data, np.log(y_data), deg=2)

型
接下来，声明一个变量X，并根据这个变量创建一个多项式。将它提升到e的幂。

x = sm.Symbol('x')
logged_poly = sm.Poly(poly_coeffs, x).as_expr()
non_logged_poly = sm.exp(logged_poly)

型
接下来，求这个表达式对x的二阶导数。

twice_diff = sm.diff(non_logged_poly, x, 2)

型
这给了我们一个表达式
x1c 0d1x的数据
这是相当丑陋的，但至少避免了使用np.gradient()来寻找二阶导数。
现在，要在代码中使用它，您需要将其转换回Python函数。SymPy函数lambdify()可以用来完成此任务。

f = sm.lambdify(x, twice_diff, 'numpy')

型
如果你想看到生成函数的源代码，可以使用help(f)来实现。（你不需要将这个源代码复制到你的程序中。你可以直接使用f()。它只是为了显示SymPy在内部做什么。）

def _lambdifygenerated(x):
    return 9.99113566549129*(0.504928487945428*(x + 0.721464680150301)**2 + 0.710583202690176)*exp(x*(0.355291601345088*x + 0.512660683049044))

型
您可以通过一次传递一个数字或传递一个数字数组来使用此函数。

f(0)
f(np.array([-1, 0, 1]))

型
如果你想绘制这个函数的结果，你可以这样做：

plt.plot(x_data, f(x_data))

型

赞(0）回复(0）举报 11个月前

我来回答

pandas Python将多项式拟合到自然对数，然后将转换反转回非对数空间？

2条答案

简单回归

日志转换

自动化

可选

求二阶导数

根据回归系数

提案

相关问题

热门标签

最新问答