回答此问题可获得 20 贡献值,回答如果被采纳可获得 50 分。
<p>我有一个数组y,它包含在每个月的某一天观察到的值。月份的日期在数组x中</p>
<p>我需要用三次样条插值这些值,这样我就可以得到一个月里每一天的值。为了考虑每月的每一天,我创建了一个数组xd。在</p>
<p>如果我想绘制原始y和插值y(即yd),我需要它们在同一个轴上对齐。这个轴是一个考虑到一个月的一天,xd。在</p>
<p>有没有一种有效的方法可以快速创建一个新的y数组,在新的x轴的正确位置精确地包含原始的y元素,所有其他元素都填充零或NaN(最好)?在</p>
<p>例如,我的第一个y只在第2天可用,所以在新的y数组中,我需要第一个元素来显示0/NaN。然后第二个元素将显示原始y=11,第三个元素将显示NaN,依此类推</p>
<p>我已经写了这段代码,上面我提到了,但我不知道是否有更好/更快的方法来实现这一点。在许多情况下,数组比我在下面的例子中展示的要大得多,所以有一些有效的算法会有所帮助。非常感谢。在</p>
<pre><code>import numpy as np
import scipy.interpolate as sp
x = [2, 5, 7, 11, 13, 16, 19, 23, 25, 30]
y = [11, 10, 12, 14, 16, 19, 17, 14, 18, 17]
xd = np.linspace(0, max(x), int(max(x))+1) # create the new x axis
ipo = sp.splrep(x, y, k=3) # cubic spline
yd = sp.splev(xd, ipo) # interpolated y values
newY = np.zeros((1, len(yd)), dtype=float) # preallocate for the filled y values
for i in x:
if(i in xd):
idx, = np.where(xd == i) # find where the original x value is in the new x axis
idx2, = np.where(np.array(x) == i)
newY[0, int(idx)] = y[int(idx2)] # replace the y value of the new vector with the y value from original set
</code></pre>
<p><strong>编辑:</strong></p>
<p>为了澄清一下,需要一组对齐的数组(两个数组共享同一个轴)是因为当我绘制两个数组(newY和yd)时,我还添加了一些子图,其中我取绝对值和相对值的差异,以查看适合性有多好。在</p>
<p>我知道在这种情况下,样条曲线将始终通过我作为输入的所有点,因此差值将为零,但下面的绘图函数应该适用于任何类型的比较(即任何类型的插值与实际输入)。我使用的绘图功能如下:</p>
^{pr2}$
<p><strong>编辑2:</strong></p>
<p>添加到目前为止提案的绩效指标。循环的速度比我的循环快得多,因为循环的速度比循环快得多。在我的例子中,x有23个元素,而xd有3655个元素。在</p>
<pre><code>def A():
for i in x:
if(i in xd):
idx, = np.where(xd == i) # find where the original x value is in the new x axis
idx2, = np.where(np.array(x) == i)
newY[int(idx)] = y[int(idx2)] # replace the y value of the new vector with the y value from original set
def B():
for i, date in enumerate(xd):
if date in x:
new_y[i] = date
def C():
known_values = dict(zip(x, y))
for i,u in enumerate(xd):
if u in known_values:
newY[i] = known_values[u]
</code></pre>
<p>%timeit A()
每个回路219µs±8.8µs(平均值±标准偏差,7次运行,每个回路1000次)</p>
<p>%timeit B()
每个回路8.87 ms±95.3µs(7次运行的平均值±标准偏差,每个回路100次)</p>
<p>%timeit C()
每个回路408µs±11.3µs(平均值±标准偏差,7次运行,每个回路1000次)</p>
<p>我还尝试将A()函数传递给Numba进行JIT编译:</p>
<pre><code>A_nb = numba.jit(A)
</code></pre>
<p>获取:</p>
<p>%timeit A_nb()
每个回路226µs±610 ns(平均值±标准偏差,7次运行,每个回路1000次)</p>