【问题标题】:scipy: interpolation, cubic & linearscipy:插值,三次和线性
【发布时间】:2023-03-04 02:42:02
【问题描述】:

我正在尝试插入我的数据集(第一列是时间,第三列是实际数据):

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d

data = np.genfromtxt("data.csv", delimiter=" ")

x = data[:, 0]
y = data[:, 2]

xx = np.linspace(x.min(), x.max(), 1000)
y_smooth = interp1d(x, y)(xx)
#y_smooth = interp1d(x, y, kind="cubic")(xx)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(xx, y_smooth, "r-")

plt.show()

但我发现线性插值和三次插值之间存在一些奇怪的差异。 这是线性的结果:

立方也是一样的:

我不确定,为什么图一直在跳,y_smooth 包含不正确的值?

ipdb> y_smooth_linear.max()
141.5481144
ipdb> y_smooth_cubic.max()
1.2663431888584225e+18

谁能给我解释一下,我怎样才能改变我的代码来实现正确的插值?

UPD:这里是data.cvs file

【问题讨论】:

标签: python scipy interpolation


【解决方案1】:

您的数据包含同一个 x 值的多个 y 值。这违反了大多数插值算法的假设。

要么丢弃具有重复 x 值的行,平均每个单独 x 的 y 值,要么为 x 值获得更好的分辨率,使它们不再相同。

【讨论】:

  • 另见scipy issuescipy.interpolate.splmake 似乎是目前scipy 中唯一支持它的插值例程。
【解决方案2】:

鉴于cfh's observation x 具有重复值,您可以使用np.unique 为每个 x 选择唯一值 y

x2, idx = np.unique(x, return_index=True)
y2 = y[idx]

return_index=True 导致np.unique 不仅返回唯一值x2,还返回原始x 数组中唯一xs 的位置idx。请注意,这会为每个唯一的 x 选择 yfirst 值。

如果您想平均每个唯一 x 的所有 y 值,您可以使用 stats.binned_statistic:

import scipy.stats as stats
x2, inv = np.unique(x, return_inverse=True)
y2, bin_edges, binnumber = stats.binned_statistic(
    x=inv, values=y, statistic='mean', bins=inv.max()+1)

return_inverse=True 告诉 np.unique 返回索引, 可以重建原始数组。这些指数也可以作为分类 标签或“因素”,这是它们在调用中的使用方式 binned_statistic 以上。


import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp1d
import scipy.stats as stats

data = np.genfromtxt("data.csv", delimiter=" ")

x = data[:, 0]
y = data[:, 1]

x2, idx, inv = np.unique(x, return_index=True, return_inverse=True)
y_uniq = y[idx]
y_ave, bin_edges, binnumber = stats.binned_statistic(
    x=inv, values=y, statistic='mean', bins=inv.max()+1)

xx = np.linspace(x.min(), x.max(), 1000)
y_smooth = interp1d(x, y)(xx)
y_smooth2 = interp1d(x2, y_uniq, kind="cubic")(xx)
y_smooth3 = interp1d(x2, y_ave, kind="cubic")(xx)

fig, ax = plt.subplots(nrows=3, sharex=True)

ax[0].plot(xx, y_smooth, "r-", label='linear')
ax[1].plot(xx, y_smooth2, "b-", label='cubic (first y)')
ax[2].plot(xx, y_smooth3, "b-", label='cubic (ave y)')
ax[0].legend(loc='best')
ax[1].legend(loc='best')
ax[2].legend(loc='best')
plt.show()

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 2018-12-27
    • 1970-01-01
    • 1970-01-01
    • 2016-06-22
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2010-12-11
    相关资源
    最近更新 更多