使用scipy.curve_fit结合mpmath.polylog拟合数据遇类型错误求助

尝试使用scipy.optimize.curve_fit结合mpmath.polylog进行数据拟合，定义拟合函数如下：

def Vendik(x, C, T0, eps_cst, Td, xsi):
    integral = (np.pi**2)/6. + (Td/x) * np.log(1-np.cosh(Td/x)+np.sinh(Td/x)) - polylog(2, np.exp(-Td/x))
    eta = ( Td/ (2*T0) ) * ( 0.5 + 2./(Td/x)**2 * integral) -1
    term1 = np.sqrt( xsi**2 + eta**3 ) + xsi
    term2 = np.sqrt( xsi**2 + eta**3 ) - xsi 
    return eps_cst + (C*1e4/T0) / ( term1**(2./3.) + term2**(2./3.) - eta )

拟合执行代码：

xdata = data1['T(K)']
ydata = data1['epsr_STO']
p1 = np.array([5.769, 13.16, 66.7, 86.5, 0.667])
popt_Vendik, pcov_Vendik = curve_fit(Vendik, xdata, ydata, p1)

触发错误：

File /opt/anaconda3/lib/python3.9/site-packages/mpmath/ctx_mp.py:634 in _convert_fallback
    raise TypeError("cannot create mpf from " + repr(x))

但通过循环逐个调用函数（res= [Vendik(data1['T(K)'][i], 5.769, 13.16, 66.7, 86.5, 0.667) for i in range(len(data1))]）无报错，且使用不依赖mpmath的拟合函数处理相同数据可正常完成。已尝试将xdata/ydata转换为np.array()或list，问题依旧。

mpmath.polylog仅支持标量输入，而curve_fit在拟合过程中会批量传递numpy数组给拟合函数。循环调用时传入的是单个标量值，因此可正常运行；但curve_fit传入数组时，mpmath.polylog无法处理数组类型，导致类型转换错误。

方案1：向量化mpmath.polylog #

用np.vectorize包装mpmath.polylog，使其支持numpy数组输入：

import numpy as np
import mpmath
from scipy.optimize import curve_fit

# 向量化处理polylog，使其支持数组
vectorized_polylog = np.vectorize(lambda n, z: mpmath.polylog(n, z))

def Vendik(x, C, T0, eps_cst, Td, xsi):
    exp_term = np.exp(-Td/x)
    integral = (np.pi**2)/6. + (Td/x) * np.log(1-np.cosh(Td/x)+np.sinh(Td/x)) - vectorized_polylog(2, exp_term)
    eta = ( Td/ (2*T0) ) * ( 0.5 + 2./(Td/x)**2 * integral) -1
    term1 = np.sqrt( xsi**2 + eta**3 ) + xsi
    term2 = np.sqrt( xsi**2 + eta**3 ) - xsi 
    return eps_cst + (C*1e4/T0) / ( term1**(2./3.) + term2**(2./3.) - eta )

方案2：改用scipy.special.polylog（推荐） #

若你的scipy版本≥0.19.0，scipy.special.polylog原生支持numpy数组，无需额外处理，直接替换即可：

import numpy as np
from scipy.special import polylog
from scipy.optimize import curve_fit

def Vendik(x, C, T0, eps_cst, Td, xsi):
    integral = (np.pi**2)/6. + (Td/x) * np.log(1-np.cosh(Td/x)+np.sinh(Td/x)) - polylog(2, np.exp(-Td/x))
    eta = ( Td/ (2*T0) ) * ( 0.5 + 2./(Td/x)**2 * integral) -1
    term1 = np.sqrt( xsi**2 + eta**3 ) + xsi
    term2 = np.sqrt( xsi**2 + eta**3 ) - xsi 
    return eps_cst + (C*1e4/T0) / ( term1**(2./3.) + term2**(2./3.) - eta )

验证 #

修改后重新执行拟合代码即可，两种方案均能解决数组输入的兼容问题，方案2性能更优（原生数组运算，无需循环处理单个元素）。

内容的提问来源于stack exchange，提问作者Jep