使用 Python 为每个数据点拟合具有多个输入、多个输出、多个参数和协方差矩阵的模型

Question

注意：类似的问题已在交叉验证站点针对数学/理论部分提出，按照社区的建议。这里的问题主要是关于Python中拟合过程的实现。

我想拟合一个具有

的模型

1. 多种输入
2. 多种输出
3. 多个参数

让我做一个简单的示例以使其更容易。例如，假设我有一个 2D 矢量场模型，我想用它来拟合一些数据。因此，该向量场是一个具有 2 个输入（

x

、

y

位置）和 2 个输出（这些位置处向量的

u

、

v

分量）的模型。该模型可以使用 3 个参数（

a

、

b

和

c

）进行调整。

下面的代码定义了手头的模型，并创建了一些构成我们测量集的噪声数据：

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def model(XY, a, b, c):
    x, y = XY
    u = a*x+b*y
    v = a*x*y-c
    return u, v

# TRUE MODEL FIELD
a, b, c = -1, 7, 11 # True model parameters
x, y = np.meshgrid(np.linspace(-5, 5, 10), np.linspace(-5, 5, 10))
u, v = model((x, y), a, b, c)

# PLOT MODEL
plt.style.use('dark_background')
plt.quiver(x, y, u, v, color = 'white')
plt.xlabel('x')
plt.ylabel('y')
plt.savefig('Model_Field.png', format = 'png', dpi = 400)
plt.close('all')

# CREATE DATA
data_size = 100
np.random.seed(0)
x_data, y_data = np.random.uniform(-5, 5, data_size), np.random.uniform(-5, 5, data_size)
u_data_expected, v_data_expected = model((x_data, y_data), a, b, c)
noise_field = [np.random.normal(0, 5, data_size), np.random.normal(0, 5, data_size)]
u_data_real, v_data_real = u_data_expected + noise_field[0], v_data_expected + noise_field[1]

# PLOT DATA
plt.quiver(x_data, y_data, u_data_real, v_data_real, color = 'crimson')
plt.quiver(x_data, y_data, u_data_expected, v_data_expected, color = 'aqua', alpha = 0.6)
plt.xlabel('x')
plt.ylabel('y')
plt.savefig('Data_Field.png', format = 'png', dpi = 400)
plt.close('all')

我们可以通过下图了解

a = -1

、

b = 7

和

c = 11

的场模型是什么样的：

第二张图显示了我们将在其上拟合模型的虚构 (

x

,

y

,

u

,

v

) 数据。蓝色是“测量”位置处的预期矢量场，红色是“测量”数据（我添加了一些噪声）：

因此，为了简单起见，我想对这些数据（上图中的红色向量）进行拟合并获得最佳拟合参数及其协方差（这对我来说很重要）。我不介意使用 Scipy、lmfit 或任何其他可能适合此任务的库。

如果可以做到这一点，考虑到我的数据具有以下事实，我还想进行拟合：

4. 输出的不确定性 (

   u_err

   ,

   v_err

   )
5. 输入的不确定性 (

   x_err

   ,

   y_err

   )
6. 输入数据之间的相关性 (

   x_y_corr

   )
7. 输出数据之间的相关性 (

   u_v_corr

   )

简而言之，每个输入坐标都有一个协方差矩阵，每个输出/测量向量都有一个协方差矩阵，在拟合模型参数时需要考虑所有这些。

这里我给出了一些数据协方差矩阵的示例：

# CREATE DATA UNCERTAINTIES AND CORRELATIONS
# Input data (positions) uncertainties and correlations:
x_err, y_err = np.random.normal(0, 0.1, data_size), np.random.normal(0, 0.1, data_size)
x_y_corr = np.random.uniform(-1, 1, data_size)
cov_matrices_xy = np.array([[[x_err[i]**2, x_err[i]*y_err[i]*x_y_corr[i]],
                             [x_err[i]*y_err[i]*x_y_corr[i], y_err[i]**2]] for i in range(data_size)])

# Output data (vectors) uncertainties and correlations:
u_err, v_err = np.random.normal(0, 0.1, data_size), np.random.normal(0, 0.1, data_size)
u_v_corr = np.random.uniform(-1, 1, data_size)
cov_matrices_uv = np.array([[[u_err[i]**2, u_err[i]*v_err[i]*u_v_corr[i]],
                             [u_err[i]*v_err[i]*u_v_corr[i], v_err[i]**2]] for i in range(data_size)])

所以，简而言之，我的问题是如何在 Python 中拟合模型并考虑列出的 1-7 个属性？

---

我尝试了什么？

嗯，我已经能够满足第 1、3 和 4 点（因此，用所有输入变量拟合模型，但仅在其中一个输出变量上拟合模型，并且在考虑该输出变量的不确定性的同时，一切都进行拟合）。

from scipy.optimize import curve_fit

# Redefine the model function so that it only yields one of the vector components:
def model(XY, a, b, c):
    x, y = XY
    u = a*x+b*y
    v = a*x*y-c
    return u

# FITTING THE MODEL
init_parameters = [1, 3, 5] # a, b, c
params_fit, cov_params_fit = curve_fit(f = model,
                                       xdata = (x_data, y_data),
                                       ydata = u_data_real,
                                       sigma = u_err,
                                       p0 = init_parameters)

得到的参数是

a = -0.78162045

（非常接近预期的

a = -1

）、

b = 6.41216698

（非常接近预期的

b = 7

）和

c = 5

（与预期的

c = 11

无关，因为拟合仅在

u

组件上进行，并且

c

参数对此不起作用，因此它保持最初的状态）。由于某种原因，我无法获得参数的协方差矩阵，

cov_params_fit

，但这可能与这个问题无关。

所以，我仍然不知道如何考虑第 2、5、6 和 7 点（因此，考虑输入的不确定性以及输入变量与输出变量之间的相关性，并且在拟合 2D 向量时数据不仅仅是组件之一）。

为了解决第 2 点，我可以将模型函数拆分为两个模型函数（每个组件一个）并分别拟合每个模型函数，但这样我就会有两组

a,b,c

参数，并且拟合效果不会像它那样好在两个组件上同时完成。

我也尝试用这段代码解决2和4：

# Redefine the model function so that it yields both components again:
def model(XY, a, b, c):
    x, y = XY
    u = a*x+b*y
    v = a*x*y-c
    return u, v

# FITTING THE MODEL
init_parameters = [1, 3, 5] # a, b, c
params_fit, cov_params_fit = curve_fit(f = model,
                                       xdata = (x_data, y_data),
                                       ydata = (u_data_real, v_data_real),
                                       sigma = (u_err, v_err),
                                       p0 = init_parameters)

但这会产生错误。看来

curve_fit()

可能无法处理如此高级的配件，但也许有人知道如何做到这一点。

解决 2（模型事物的多个输出）的另一种想法可以将模型函数转换为仅返回一个复值输出的函数，

u+vj

，并执行类似于数据向量的操作。也许数据可以与仅适用于一个输出值的程序相匹配。

Answer 1

只需修改模型的签名即可堆叠字段维度。这可以解决你的问题。

def fitter(XY, a, b, c):
    x, y = XY
    u = a * x + b * y
    v = a * x * y - c
    return np.hstack([np.ravel(u), np.ravel(v)])

然后你可以回归你的数据集：

target = np.hstack([np.ravel(u_data_real), np.ravel(v_data_real)])
optimize.curve_fit(fitter, (x_data, y_data), target)

# (array([-1.05860421,  6.86975359, 11.70709473]),
#  array([[ 3.45604885e-03,  2.69747434e-04, -2.08667038e-03],
#         [ 2.69747434e-04,  3.03704636e-02, -1.62866327e-04],
#         [-2.08667038e-03, -1.62866327e-04,  2.36032115e-01]]))

这似乎与你的模型一致，但有错误。