返回numpy数组的函数的曲线拟合

Question

我知道{}的{}库及其拟合曲线的能力。我在这里和文档中读了很多例子，但我无法解决我的问题。 例如，我有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, epsilon, 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中描述的方法那样解决此问题

Answer 1

曲线拟合

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

[It uses] non-linear least squares to fit a function, f, to data.

它指的是least_squares，在这里它表示

Given the residuals f(x) (an m-D real function of n real variables) and the loss function rho(s) (a scalar function), least_squares finds a local minimum of the cost function 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`

凸性

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

1. Lennard-Jones函数具有排斥力和吸引力之差，两者在距离上均为负偶数指数，仅此一项不太可能是凸的

2. 以不同位置为中心的多个局部函数之和没有定义的凸性

3. 众所周知，分子能、晶体能或蛋白质折叠是非凸的

几天前（骑自行车时）我在想这个问题，分子如何在全球最小能量下配置，我想知道它是否因为量子隧穿效应而发现配置如此迅速

非凸优化

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

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

您将在scipy.optimize中找到一些global optimization方法。我鼓励您尝试Basin hopping，因为它成功地应用于类似的问题，正如您在references中看到的那样

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

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