python - 如何使用 numpy 在 2D 网格上矢量化卡方计算？

Question

N我有一些由 3 个长度为: x, y,的一维数组组成的数据y_errors，我正在使用A*sin(w*x+phi)+c. 我已经定义了一个卡方函数并用于scipy.optimise.minimise将模型拟合到数据中。我现在想考虑一个离散的二维网格，范围覆盖来自A, w,的任何 2 个参数的一组值phi，c例如，我可能想考虑A范围超过[0,1,2,3]和w范围超过的值的网格[0.2,0.3,0.4]。我想计算网格上每对值的卡方值（使用网格未跨越的参数的固定值，即phi在c示例中）。

这可以使用 for 循环轻松实现，但效率很低。我想找到一种更快的方法来进行计算。我对其他想法持开放态度，但最明显的解决方案似乎是使用 numpy 数组对问题进行矢量化。我的问题是：我如何最有效地（最好是直接地）计算网格上的所有卡方值（使用 numpy 向量化方程）？

在堆栈交换中也有类似的问题，例如Python numpy 中的无环卡方网格搜索以及如何对这个特定的非 numpy 函数进行矢量化？，但是对于给出最佳/最有效或“正确”的方法，答案似乎并不十分自信。我试图解决的问题是物理科学中数据分析的一个非常普遍和重要的任务（仅我大学就有数百名本科生会喜欢一个解决方案），所以我认为值得在上下文中得到一个明确而明确的答案我认为这个例子是解决这类问题的一个很好的典型设置。

我试图自己解决这个问题，但是在将数组投射在一起时遇到了困难。我的方法是为可变参数（例如Aand w）的 2D 网格制作一个 numpy meshgrid，然后我将这些与常量参数（例如phiand c）一起存储在一个列表中，我将其变成一个 numpy 数组和这个参数数组提供给卡方，因此提供给模型函数以进行卡方计算（我还尝试制作与网格网格中的每个数组相同形状的 numpy 数组，并填充剩余 2 个参数的常量值 - 这在我的代码中被注释掉了）。我遇到的问题是，在计算卡方时，这需要遍历数组中的所有数据点x, y,y_errors以及迭代网格网格中的所有不同参数值。尝试这样做会导致 numpy 转换错误。

我将在下面包含我的代码的相关部分以供参考：

import pandas as pd
import numpy as np
import scipy.optimize as opt
import scipy.stats as stats

#TEST DATA
x=np.linspace(0,2*np.pi,20)
y=np.sin(x)+np.random.normal(loc=0,scale=0.01,size=len(x))
y_error=np.array([0.01]*len(x))

#MODEL
def sinusoidal_model(t,params):
    A,w,phi,c=params
    return(A*np.sin(w*t+phi)+c)

#CHI-SQUARED
def chi_squared(model_params, model, x, y, y_error):
    return np.sum(((y - model(x, model_params))/y_error)**2)

#INITIAL MODEL PARAMETERS
A=1
w=1
phi=0
c=0
initial_params=np.array((A,w,phi,c))

#FITTING
model=sinusoidal_model
deg_freedom = x.size - initial_params.size 
fit = opt.minimize(chi_squared, initial_params, args=(model, x, y, y_error))
fit_params= fit.x
chisq_min = chi_squared(fit_params, model, x, y, y_error)

#GENERATE COLOUR AND CONTOUR PLOTS
i_of_param_i=0
j_of_param_j=1
param_i=fit_params[i_of_param_i]
param_j=fit_params[j_of_param_j]
param_i_half_range=1
param_j_half_range=1
n_points_x=5
n_points_y=4
param_i_axis=np.linspace(param_i-param_i_half_range,param_i+param_i_half_range,n_points_x)
param_j_axis=np.linspace(param_j-param_j_half_range,param_j+param_j_half_range,n_points_y)
X,Y=np.meshgrid(param_i_axis,param_j_axis,indexing='xy')
"""
grid_shape=X.shape
parameters=[np.full(grid_shape,fit_params[i]) for i in range(len(fit_params))]
parameters[i_of_param_i]=X
parameters[j_of_param_j]=Y
parameters=np.array(parameters)
x=np.full(grid_shape,x)
y=np.full(grid_shape,y)
y_error=np.full(grid_shape,y_error)
Z=chi_squared(parameters, model, x, y, y_error)
"""
parameters=list(np.array(param) for param in list(fit_params))
parameters[i_of_param_i]=X
parameters[j_of_param_j]=Y
parameters=np.array(parameters)
Z=chi_squared(parameters, model, x, y, y_error)

请注意，“生成颜色和轮廓图”部分的前几行就在那里，以便我可以轻松地切换我想要保持不变的参数以及我想要在网格上改变的参数。

摘要：我想在一些 2D 参数空间上为模型制作卡方曲面图，而其他任何参数都保持固定。问题是卡方必须考虑数据和计算卡方的参数范围，但这些有不同的维度，不能相互转换。

标签： pythonnumpyvectorizationchi-squared