为什么我在 scipy curve_fit 中的近似值如此糟糕？

Question

在对我的计算数据进行近似计算之前，我决定在 Q 的数据上测试代码。 Zhang, N. Sabelli, V. Buch； H···H2O 的势能面。 J.化学。物理。 1991 年 7 月 15 日； 95 (2): 1080–1085，公共领域的一篇旧文章，遇到了一个问题：我用文章中相同的函数进行的逼近结果与数据和文章中的逼近都相差甚远，它甚至没有最小值，这对于这一系列任务非常重要

在画文章近似图的时候，我舍弃了每个几何的第一个点（除了最后一个，文章近似和我的都不好），因为能量很大（显然是为了在中有一个很好的对应关系）最小面积）和图表结果是无用的。

云与data.

from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import math
import scipy
from numpy import array

data = pd.read_excel('data.xlsx', header=None)
values_E = array(data[3])
values_R = array(data[0])
values_theta = array([math.radians(data[1][i]) for i in range(len(data[1]))])
values_phi = array([math.radians(data[2][i]) for i in range(len(data[2]))])
values = array([values_R, values_theta, values_phi])

def mapping(values, eps00, eps10, eps20, eps2_2, sigma00, sigma10, sigma20, sigma2_2):
    return [math.fsum([
    4*eps00*(sigma00/values[0][i])**6*((sigma00/values[0][i])**6 - 1)* 1/2 * np.sqrt(1/np.pi),
    4*eps10*(sigma10/values[0][i])**6*((sigma10/values[0][i])**6 - 1)* 1/2 * np.sqrt(3/np.pi) * np.cos(values[1][i]),
    4*eps20*(sigma20/values[0][i])**6*((sigma20/values[0][i])**6 - 1)* 1/4 * np.sqrt(5/np.pi) * (3 * (np.cos(values[1][i]))**2 - 1),
    4*eps2_2*(sigma2_2/values[0][i])**6*((sigma2_2/values[0][i])**6 - 1)* 1/4 * np.sqrt(15/np.pi) * (np.sin(values[1][i]))**2 * (np.sin(values[2][i]))**2
    ]) for i in range(len(values_R))]
param_bounds=([0,-20, 0, 0, 0, 0, 0, 0], [150,10,10,100,10,10,10,10])
args, _ = curve_fit(mapping, values, values_E, bounds=param_bounds, maxfev=10000000)
eps00, eps10, eps20, eps2_2, sigma00, sigma10, sigma20, sigma2_2 = args

#energy by approximation from the article and by my
article = mapping(values, 107.9, -9.5, 8.6, 35.6, 3.0, 3.3, 2.98, 2.92) 
new_E = mapping(values, eps00, eps10, eps20, eps2_2, sigma00, sigma10, sigma20, sigma2_2)

sp = plt.subplot(231) #90_00
plt.plot(values_R[0:14], values_E[0:14], label='E_90_00')
plt.plot(values_R[1:14], article[1:14], label='article')
plt.plot(values_R[0:14], new_E[0:14], label='new_E')
plt.xlabel('R_90_00')
plt.ylabel('E_90_00')
plt.legend(loc='best')

sp = plt.subplot(232) #90_90
plt.plot(values_R[14:28], values_E[14:28], label='E_90_90')
plt.plot(values_R[15:28], article[15:28], label='article')
plt.plot(values_R[14:28], new_E[14:28], label='new_E')
plt.xlabel('R_90_90')
plt.ylabel('E_90_90')
plt.legend(loc='best')

sp = plt.subplot(233) #00_00
plt.plot(values_R[28:42], values_E[28:42], label='E_00_00')
plt.plot(values_R[29:42], article[29:42], label='article')
plt.plot(values_R[28:42], new_E[28:42], label='new_E')
plt.xlabel('R_00_00')
plt.ylabel('E_00_00')
plt.legend(loc='best')

sp = plt.subplot(234) #180_00
plt.plot(values_R[42:56], values_E[42:56], label='E_180_00')
plt.plot(values_R[43:56], article[43:56], label='article')
plt.plot(values_R[42:56], new_E[42:56], label='new_E')
plt.xlabel('R_180_00')
plt.ylabel('E_180_00')
plt.legend(loc='best')

sp = plt.subplot(235) #127_90
plt.plot(values_R[56:62], values_E[56:62], label='E_127_90')
plt.plot(values_R[56:62], article[56:62], label='article')
plt.plot(values_R[56:62], new_E[56:62], label='new_E')
plt.xlabel('R_127_90')
plt.ylabel('E_127_90')
plt.legend(loc='best')

plt.show()

我检查了输入并重写了映射以更好看，并对参数做了非常精确的限制，当你不知道结果是什么时这是不可能的，但这一切都是徒劳的。 我没有描述原始问题和文章，因为目前一切都归结为代码问题。

我的能量曲线应该与文章中的数据曲线重合，或者至少与文章中的近似值重合。

Answer 1

大量不同严重程度的错误：

• 当你从这个电子表格中读入 Pandas 数据框时，你需要提供列名；不要使用位置列索引
• 不要在任何地方使用

  math

在考虑文章的术语时，

• mapping

  不是您函数的合适名称

• fsum

  是关心准确性的错误地方，应该离开
• 向量化的普遍失败
• 未能在您的绘图中正确使用数据框分组

但是所有这些仍然不能解释您的大部分错误。我有根据的猜测如下：

请注意，在表 III 中的文章中，所有角度组的半径从 2 到 6 埃，但最后一组 127.74° / 90° 除外。我相信他们对那个群体的适应性比你好得多，而他们对所有其他群体的适应性都比你差，这绝非巧合。我认为发生的事情是，他们的拟合应用于比他们在文章中打印的更大的角度组，或，他们应用了与您不同的拟合误差权重，以便它们根据每个组的幅度进行归一化。

这个理论得到了以下事实的证实：如果你只适合这个组，你会看到与他们相似的效果：那个组的适合度变得更好，而其他每个组的适合度都有一个太高的起始峰值。

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

data = pd.read_csv('data.csv', header=None, names=('R', 'theta', 'phi', 'E'), decimal=',')
data[['theta_rad', 'phi_rad']] = np.deg2rad(data[['theta', 'phi']])
values = data[['R', 'theta_rad', 'phi_rad']].values.T


def Vfit(
    values: np.ndarray,
      eps0p0: float,   eps1p0: float,   eps2p0: float,   eps2n2: float,
    sigma0p0: float, sigma1p0: float, sigma2p0: float, sigma2n2: float,
) -> np.ndarray:
    R, theta, phi = values

    cos_theta = np.cos(theta)
    Y1p0 = 0.50*np.sqrt( 3/np.pi) * cos_theta
    Y2p0 = 0.25*np.sqrt( 5/np.pi) * (3*cos_theta**2 - 1)
    Y2n2 = 0.25*np.sqrt(15/np.pi) * np.sin(theta)**2 * np.sin(phi)**2
    Y0p0 = np.full_like(Y1p0, fill_value=0.50*np.sqrt(1/np.pi))  # 1*pi?

    Y = np.array((Y0p0, Y1p0, Y2p0, Y2n2))
    eps = (eps0p0,), (eps1p0,), (eps2p0,), (eps2n2,),
    sigma = (sigma0p0,), (sigma1p0,), (sigma2p0,), (sigma2n2,),

    ratio = sigma/R
    product = 4*Y*eps*(ratio**12 - ratio**6)
    return product.sum(axis=0)


is_fit = (data.phi.astype(int) == 90) & (data.theta.astype(int) == 127)

args, _ = curve_fit(
    f=Vfit, xdata=values[:, is_fit], ydata=data.E[is_fit],
    bounds=(
        # epsilon        sigma
        (  0,-20, 0,  0,  0, 0, 0, 0),
        (150, 10,10,100, 10,10,10,10),
    ),
    p0=(110, -10, 10, 30, 3, 3, 3, 3),
    maxfev=10_000,
)
eps = args[:4]
sigma = args[4:]
print('Fit epsilon:', eps)
print('Fit sigma:', sigma)

params_fbd = 107.9, -9.5, 8.6, 35.6, 3.0, 3.3, 2.98, 2.92
article = Vfit(values, *params_fbd)
fit_E = Vfit(values, *eps, *sigma)

fig, ax_rows = plt.subplots(nrows=2, ncols=3)
all_axes = [ax for row in ax_rows for ax in row]
for ax, ((theta, phi), group) in zip(all_axes, data.groupby(['theta', 'phi'])):
    ax.plot(group.R, group.E, label='E')
    ax.plot(group.R, article[group.index], label='article')
    ax.plot(group.R, fit_E[group.index], label='fit_E')
    ax.set_xlabel('R (Å)')
    ax.set_ylabel('E (cm⁻¹)')
    ax.set_title(f'theta={theta} phi={phi}')
    ax.legend()

plt.show()

Fit epsilon: [113.92910577 -11.10362669   9.99792662  38.44309601]
Fit sigma: [3.09518901 3.09459379 3.09304591 3.09410421]