python - 用面具Python去趋势

我有一个file我读到如下：

data1 = np.loadtxt('lc1.out')
x = data1[:, 0]
y = data1[:, 1]

我想删除它，我发现了一个非常有用的链接here。

model = np.polyfit(x, y, 2)
predicted = np.polyval(model, x)

无论如何，我想屏蔽一部分数据，例如我将只使用屏蔽外的点来拟合。例如，我只想使用小于639.5且大于641.5的数据和二阶多项式拟合。
我有使用ma.masked_outside(x, 639.5, 641.5)的想法，因为它将很容易保存在一个数组中，只保存掩码之外的元素…但我不知道如何使用polyfit来转换它。

最佳答案

在您的用例中使用屏蔽数组可能没有什么困难的原因，除非是出于性能原因或是为了进一步使用屏蔽。所以我将演示如何使用和不使用蒙面阵列。
但我们先不加掩饰地进行描述，以便获得参考：
无掩蔽的二阶多项式分解

import numpy as np
import matplotlib.pyplot as plt

data1 = np.loadtxt('lc1.out')
x, y = data1.T

fig = plt.figure()

plt.subplot(2, 1, 1)
plt.title('polyfit, original data set')
plt.plot(x, y, 'c.')

coeff = np.polyfit(x, y, 2)

# no need to use the original x values here just for visualizing the polynomial
x_poly = np.linspace(x.min(), x.max())
y_poly = np.polyval(coeff, x_poly)
plt.plot(x_poly, y_poly, 'r-', linewidth=3)

mid = len(x_poly) // 2
plt.annotate('y = {:.7g} x\xB2 + {:.7g} x + {:.7g}'.format(*coeff),
             (x_poly[mid], y_poly[mid]), (0, 48), textcoords='offset points',
             arrowprops={'arrowstyle': '->'}, horizontalalignment='center')

plt.subplot(2, 1, 2)
plt.title('detrended')

# we need the original x values here, so we can remove the trend from all points
trend = np.polyval(coeff, x)
# note that simply subtracting the trend might not be enough for other data sets
plt.plot(x, y - trend, 'b.')
fig.show()

注意多项式的系数。
二阶多项式分解，选择有效点
我们可以简单地创建新的x和y数组，其中只包含所需的点。这里出错的可能性较小。
这有三个步骤。首先，我们在感兴趣的数组上使用比较运算符，这将导致bool数组的索引处有“True”值，而其他地方的索引处都有“False”值。
然后我们将bool数组放到'np.where（）'中，结果得到一个包含所有索引号的数组，作为bool数组具有'True'值的值，也就是说，它回答了一个问题：“我的数组truthy在哪里？”
最后，我们仔细研究了NUMPY的高级索引，并将索引的结果数组应用到X和Y数组中，从而过滤掉所有不需要的索引。

import numpy as np
import matplotlib.pyplot as plt

data1 = np.loadtxt('lc1.out')
x, y = data1.T
select = np.where((x < 640.75) | (x > 641.25))
x_selection = x[select]  # numpy advanced indexing
y_selection = y[select]  # numpy advanced indexing

fig = plt.figure()

plt.subplot(2, 1, 1)
plt.title('polyfit, selecting significant points')
plt.plot(x_selection, y_selection, 'c.')

coeff = np.polyfit(x_selection, y_selection, 2)

# no need to use the original x values here just for visualizing the polynomial
x_poly = np.linspace(x_selection.min(), x_selection.max())
y_poly = np.polyval(coeff, x_poly)
plt.plot(x_poly, y_poly, 'r-', linewidth=3)

mid = len(x_poly) // 2
plt.annotate('y = {:.7g} x\xB2 + {:.7g} x + {:.7g}'.format(*coeff),
             (x_poly[mid], y_poly[mid]), (0, 48), textcoords='offset points',
             arrowprops={'arrowstyle': '->'}, horizontalalignment='center')

plt.subplot(2, 1, 2)
plt.title('detrended')

# we need the original x values here, so we can remove the trend from all points
trend = np.polyval(coeff, x)
# note that simply subtracting the trend might not be enough for other data sets
plt.plot(x, y - trend, 'b.')
fig.show()

正如所料，系数现在不同了。
二阶多项式去趋势，屏蔽不需要的点
当然，我们也可以使用蒙面数组。注意相反的逻辑：蒙面点是我们不想要的。在示例数据中，我们不需要区间内的点，我们使用ma.masked_inside()。
如果出于性能原因，我们希望避免创建原始数组的完整副本，则可以使用关键字copy=False。将原始数组设置为只读可以避免通过更改原始数组而意外更改屏蔽数组中的值。
对于屏蔽数组，我们需要使用polyfit()子模块中numpy.ma函数的版本，该版本正确地忽略了x的屏蔽版本和y的无屏蔽伴生数组中不需要的值。如果这样做失败，我们将得到错误的结果。请注意，这是一个容易犯的错误，所以如果我们可以帮助的话，我们可能会坚持使用另一种方法。

import numpy as np
import numpy.ma as ma
import matplotlib.pyplot as plt

data1 = np.loadtxt('lc1.out')
x, y = data1.T
x.flags.writeable = False  # safety measure, as we don't copy
x_masked = ma.masked_inside(x, 640.75, 641.25, copy=False)

fig = plt.figure()

plt.subplot(2, 1, 1)
plt.title('polyfit, masking unwanted points')
plt.plot(x_masked, y, 'c.')

coeff = ma.polyfit(x_masked, y, 2)

# no need to use the original x values here just for visualizing the polynomial
x_poly = np.linspace(x_masked.min(), x_masked.max())
y_poly = np.polyval(coeff, x_poly)
plt.plot(x_poly, y_poly, 'r-', linewidth=3)

mid = len(x_poly) // 2
plt.annotate('y = {:.7g} x\xB2 + {:.7g} x + {:.7g}'.format(*coeff),
             (x_poly[mid], y_poly[mid]), (0, 48), textcoords='offset points',
             arrowprops={'arrowstyle': '->'}, horizontalalignment='center')

plt.subplot(2, 1, 2)
plt.title('detrended')

# we need the original x values here, so we can remove the trend from all points
trend = np.polyval(coeff, x)
# note that simply subtracting the trend might not be enough for other data sets
plt.plot(x, y - trend, 'b.')
fig.show()

系数与其他方法相同，效果良好。如果我们错误地使用了np.polyfit()，那么我们最终得到的系数将与未屏蔽引用中的系数相同。

关于python - 用面具Python去趋势，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/43522475/