如何在部分数据点为NaN的情况下拟合二维函数？

Question

我想要将数据拟合成二维表面，具体来说，我希望找到一个将像素坐标映射到波长坐标的函数，就像IRAF中的FITCOORDS一样。

例如，在以下代码片段中，我想要找到对test数组进行拟合的结果：

import numpy as np
from astropy.modeling.models import Chebyshev2D
from astropy.modeling.fitting import LevMarLSQFitter
#%%
test = np.array([[7473, 7040, 6613, 6183, 5753, 5321, 4888],
   [7474, 7042, 6616, 6186, np.nan, 5325, 4893],
   [7476, 7044, 6619, 6189, 5759, 5328, 4897],
   [7479, 7047, np.nan, 6192, 5762, 5331, 4900]])
grid_pix, grid_wave = np.mgrid[:4, :7]
fitter = LevMarLSQFitter()
c2_init = Chebyshev2D(x_degree=3, y_degree=3)
c2_fit = fitter(c2_init, grid_wave, grid_pix, test)
print(c2_fit)

在Python 3.6上，使用astropy 2.0.2和numpy 1.13.3的结果如下：

Model: Chebyshev2D
Inputs: ('x', 'y')
Outputs: ('z',)
Model set size: 1
X-Degree: 3
Y-Degree: 3
Parameters:
    c0_0 c1_0 c2_0 c3_0 c0_1 c1_1 c2_1 c3_1 c0_2 c1_2 c2_2 c3_2 c0_3 c1_3 c2_3 c3_3
    ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
     0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
WARNING: Model is linear in parameters; consider using linear fitting methods. [astropy.modeling.fitting]

很明显，拟合从未成功过。

如果我将np.nan更改为某些值，则拟合按预期工作（例如，手动将np.nan更改为0、1等）。

我应该如何获得合理的结果？我可以让拟合器忽略np.nan值吗？

- ysBach

---

只是删除nan值而不是用虚拟值替换它们？ - Julien

{btsdaf} - ysBach

np.where 例如 - Julien

{btsdaf} - ysBach

Answer 1

您可以简单地删除数据中的nan和网格中相应的“索引”。例如，使用布尔索引：

notnans = np.isfinite(test)  # array containing True for finite values and False for nans/infs
c2_fit = fitter(c2_init, grid_wave[notnans], grid_pix[notnans], test[notnans])
print(c2_fit)

它仍会打印有关参数中的线性警告，但值确实为非零：

Model: Chebyshev2D
Inputs: ('x', 'y')
Outputs: ('z',)
Model set size: 1
X-Degree: 3
Y-Degree: 3
Parameters:
         c0_0          c1_0           c2_0      ...       c2_3            c3_3      
    ------------- -------------- -------------- ... --------------- ----------------
    7473.01546325 -431.633443323 0.471726190475 ... 0.0229037267082 -0.0012077294686

这里的诀窍在于，x、y和你的data不需要是2D数组，只要它们“表示”一个2D网格就可以是1D数组（由布尔索引返回）。

如果你有包含NaN值的“大区域”，这种方法可能不够好，因为拟合器可以适应任何东西。如果是这种情况，你可以使用astropy.convolution.convolve进行插值，然后用convolve的结果替换你的data中的NaN值。

from astropy.convolution import convolve
# Just for illustration I used a 5x5 mean filter here, the size must be adjusted
# depending on the size of all-nan-regions in your data.
mean_convolved = convolve(test, np.ones((5, 5)), boundary='extend')
# Replacing NaNs in test with the mean_convolved values:
nan_mask = np.isnan(test)
test[nan_mask] = mean_convolved[nan_mask]

# Now pass test to the fitter:
c2_fit = fitter(c2_init, grid_wave, grid_pix, test)

然而，对于一些稀疏分布的NaN值，不需要进行卷积操作。您可能需要比较这两种方法，并查看哪种方法提供更可靠的结果。缺失值可能会在拟合时造成实际问题。

- MSeifert

---

哦！这个太简单了，而且正是我想要的！非常感谢你:D - ysBach