Curve_fit 用于返回 numpy 数组的函数

如何解决Curve_fit 用于返回 numpy 数组的函数

我知道 curve_fit 的库 scipy 及其拟合曲线的能力。我在这里和文档中阅读了许多示例，但我无法解决我的问题。例如，我有 10 个文件（化学结构，但无所谓）和 10 个实验能量值。我在一个类中有一个函数，它为每个结构计算某些参数的理论能量，并返回一个包含理论能量值的 numpy 数组。

我想找到理论值最接近实验值的最佳参数。我将在这里提供我的代码的最小示例

这是读取实验能量文件、提取正确子串并将值作为numpy数组返回的类函数。 self.path 只是目录和 self.nPoints = 10。这不是那么重要，但为了完整起见我提供

def experimentalValues(self):
        os.chdir(self.path)
        energy = np.zeros(self.nPoints)
        for i in range(1,self.nPoints):
            f = open("p_" + str(i + 1) + ".xyz","r")
            energy[i] = float(f.readlines()[1].split()[1])
            f.close()
        os.chdir('..')
        return energy

我用这个以两个 numpy 数组作为参数的类函数计算了理论值，比如说

sigma = np.full(nSubstrate,2.)
epsilon = np.full(nSubstrate,0.15)

哪里nSubstrate = 9

这里是类函数。它读取文件并执行两个嵌套循环来为每个文件计算理论值并将其返回到一个 numpy 数组。

def theoreticalEnergy(self,epsilon,sigma):
        os.chdir(self.path)
        cE = np.zeros(self.nPoints)
        for n in range(0,self.nPoints):
            filenameXYZ = "p_" + str(n + 1) + "_extended.xyz"

            allCoordinates = np.loadtxt(filenameXYZ,skiprows = 0,usecols = (1,2,3))
            substrate = allCoordinates[0:self.nSubstrate]
            surface = allCoordinates[self.nSubstrate:]
            for i in range(0,substrate.shape[0]):
                positionAtomI = np.array(substrate[i][:])
                for j in range(0,surface.shape[0]):
                    positionAtomJ = np.array(surface[j][:])
                    distanceIJ = self.distance(positionAtomI,positionAtomJ)
                    cE[n] += self.LennardJones(distanceIJ,epsilon[i],sigma[i])
                
        os.chdir('..')
        return cE

同样，为了完整起见，Lennard Jones 类函数定义为

def LennardJones(self,distance,sigma):
        repulsive = (sigma/distance) ** 12.
        attractive = (sigma/distance) ** 6.
        potential = 4. * epsilon* (repulsive - attractive)
        return potential

在这种情况下，所有参数都是标量作为返回值。总结问题陈述，我有 3 个要素：

一个带有实验数据的 numpy 数组
两个 numpy 数组，猜测参数 sigma 和 epsilon
一个函数，它接受最后一个参数并返回一个包含要拟合的值的 numpy 向量。

如何像文档 https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html 中描述的方法一样解决这个问题？

解决方法

曲线拟合

curve_fit 通过找到使 f(w,x[i]) 最小化的 y[i] 使函数 w 拟合到点 sum((f(w,x[i] - y[i])**2 for i in range(n))。正如您将在函数定义后的第一行中读到的那样

[它使用] 非线性最小二乘法将函数 f 拟合到数据。

它指的是 least_squares 的地方

给定残差 f(x)（n 个实变量的 mD 实函数）和损失函数 rho(s)（标量函数），least_squares 找到成本函数的局部最小值 F(x)：

曲线拟合是一种凸成本多目标优化。由于每个单独的成本都是凸的，您可以将所有成本相加，这仍然是一个凸函数。请注意，决策变量（要优化的参数）在每个点都是相同的。

你的问题

根据我对每个能级的理解，你有一组不同的参数，如果你把它写成曲线拟合问题，目标函数可以表示为 sum((f(w[i],x[i]) - y[i])**2 ...),where y[i]is determined by the energy level. Since each of the terms in the sum is independent on the other terms,this is equivalent to finding each group of parametersw [i]separately minimizing(f(w[i],x[i]) - y[i])**2`。

凸性

凸性是一种非常方便的优化属性，因为它确保参数空间中只有一个最小值。我没有做详细的分析，但对你的能量函数的凸性有合理的怀疑。

Lennard Jones 函数具有排斥力和吸引力之差，两者在距离上都具有负偶数指数，仅此一项不太可能是凸的。
以不同位置为中心的多个局部函数的和没有定义的凸性。
众所周知，分子能量、晶体能量或蛋白质折叠是非凸的。

几天前（骑自行车）我在想这个问题，分子将如何在全局最小能量下配置，我想知道它是否会因为量子隧道效应而如此迅速地找到这种配置。

非凸优化

非凸（全局）优化与（非线性）最小二乘法不同，从某种意义上说，当找到局部最小值时，过程不会立即返回，而是开始在不同区域进行新的尝试搜索空间。如果函数是平滑的，你仍然可以利用基于梯度的局部优化方法，但复杂度仍然是 NP。

一个经典的全局优化方法是 Simulated annenaling，如果你有化学背景，我想你会对它有一些见解。 Once upon a time，在 scipy.optimize 中提供了模拟退火。

您会在 scipy.optimize 中找到一些 global optimization 方法。我鼓励您尝试 Basin hopping，因为它已成功应用于类似的问题，您可以在 references 中阅读。

我希望这能让您找到解决方案的正确方法。但是，请注意，您可能需要花时间学习如何使用该功能，并且需要做出一些决定。您需要在准确性、简单性和效率之间找到平衡点。

如果您想要更好的解决方案，请花时间推导出成本函数的梯度（您可以返回两个值 f 和 df，其中 df 是 f 相对于决策变量的梯度）。