使用curve_fit对曲线进行求和拟合

Question

我在用scipy拟合曲线时遇到了问题。这是我想要拟合我的数据的方程式：

​

​

代码如下：

import numpy as np
import os
from mpmath import *
from pandas import read_csv
from scipy.optimize import curve_fit
from matplotlib import pyplot
from sympy import *


M = 0.0129
L = 0.006
PI = np.pi

mp.dps=10


data = np.genfromtxt('Daten.csv',dtype=float,delimiter=',')

x = data[1:,0]
y = data[1:,0]

def objective(x,Mps,B,k,D):

    return (M-Mps)*(1-np.exp(-k*x))+ Mps*(-np.exp(-B*x)*np.power(D/(B*L*L),0.5)*np.power(np.tan((B*L*L)/D),0.5)- (8/(PI*PI))*\

nsum(lambda n: np.exp(-(2*n+1)*(2*n+1)*PI*PI*D*x/(4*L*L))/((2*n+1)*(2*n+1)*(1-(2*n+2)*(2*n+2)*(D*PI*PI/(4*B*L*L)))) ,[0,float('inf')]))


def obtainpar(x,y)
    popt,_ = curve_fit(objective,x,y)

    print '{Mps,B,k,D}=', popt

然后，我获得以下消息：

Traceback (most recent call last):
  File "C:/Users/gerardo.salazar/Desktop/hector.py", line 29, in <module>
    popt,_ = curve_fit(objective,x,y)
  File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\optimize\minpack.py", line 763, in curve_fit
    res = leastsq(func, p0, Dfun=jac, full_output=1, **kwargs)
  File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\optimize\minpack.py", line 388, in leastsq
    shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
  File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\optimize\minpack.py", line 26, in _check_func
    res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
  File "C:\ProgramData\Anaconda3\lib\site-packages\scipy\optimize\minpack.py", line 463, in func_wrapped
    return func(xdata, *params) - ydata
  File "C:/Users/gerardo.salazar/Desktop/hector.py", line 27, in objective
    nsum(lambda n: np.exp(-(2*n+1)*(2*n+1)*PI*PI*D*x/(4*L*L))/((2*n+1)*(2*n+1)*(1-(2*n+2)*(2*n+2)*(D*PI*PI/(4*B*L*L)))) ,[0,float('inf')]))
  File "C:\ProgramData\Anaconda3\lib\site-packages\mpmath\calculus\extrapolation.py", line 1718, in nsum
    return +ctx.adaptive_extrapolation(update, emfun, options)
  File "C:\ProgramData\Anaconda3\lib\site-packages\mpmath\calculus\extrapolation.py", line 1165, in adaptive_extrapolation
    update(partial, xrange(index, index+step))
  File "C:\ProgramData\Anaconda3\lib\site-packages\mpmath\calculus\extrapolation.py", line 1706, in update
    psum = psum + g(ctx.mpf(k))
  File "C:\ProgramData\Anaconda3\lib\site-packages\mpmath\calculus\extrapolation.py", line 1753, in g
    return f(*args)
  File "C:\ProgramData\Anaconda3\lib\site-packages\mpmath\calculus\extrapolation.py", line 1808, in g
    return f(*args)
  File "C:/Users/gerardo.salazar/Desktop/hector.py", line 27, in <lambda>
    nsum(lambda n: np.exp(-(2*n+1)*(2*n+1)*PI*PI*D*x/(4*L*L))/((2*n+1)*(2*n+1)*(1-(2*n+2)*(2*n+2)*(D*PI*PI/(4*B*L*L)))) ,[0,float('inf')]))
TypeError: loop of ufunc does not support argument 0 of type mpf which has no callable exp method

我不能分享数据，但它只是一个简单的(x，y)点，没有什么特别的。谢谢你的帮助和建议，我怎样才能做到这一点。

Answer 1

正如@piterbarg所提到的，混合使用numpy和mpmath有几个问题。我的想法是继续使用mpmath，并在计算值后切换到numpy。curve_fit等发送numpy数组，这样我就必须处理可迭代对象。最后，事实证明，负D和B是一个问题。在函数内部使用绝对值，而不是设置边界，这当然是一个选择。结果可能是，因此，negative...but不是。

相应地，代码稍微进行了清理和增强，如下所示：

import matplotlib.pyplot as plt
import numpy as np
from mpmath import pi as PI
from mpmath import nsum
from mpmath import exp
from mpmath import sqrt
from mpmath import tan
from scipy.optimize import curve_fit
from matplotlib import pyplot

M = 0.0129
L = 0.006

### might be of interest to:
# ~import mpmath
# ~mpmath.dps=10

def objective( x, Mps, b, k, d ):
    if isinstance( x, ( list, tuple, np.ndarray ) ):
        out = np.fromiter( 
            ( objective( item, Mps, b, k, d ) for item in x ), 
            np.float 
        ) ## always return np.array
    else:
        D = abs( d ) ### fails for neg
        B = abs( b ) ### fails for neg
        dil2 = D / L**2
        out = ( M - Mps ) * ( 1 - exp( -k * x ) )
        out += Mps * ( 
            -exp( -B * x ) * sqrt( dil2 / B ) 
            * sqrt( tan( B / dil2 ) ) 
            -8 / PI**2 
            * nsum( lambda n: 
                exp( -( 2 * n + 1 )**2 * PI**2 * dil2 * x / 4 ) / (
                    ( 2 * n + 1 )**2 * ( 1 - ( 2 * n + 2 )**2
                    * ( dil2 * PI**2 / ( 4 * B ) ) )
                ), [ 0, float( 'inf' ) ]
            )
        )
    return out

############ test
xl = np.linspace( -3, 15, 55 )
yl = objective( xl, 0.5, 1, 0.3, 1.2 )
yl *= np.random.normal( loc=1, scale=0.1, size=len( yl ) )


popt, pcov = curve_fit( objective, xl, yl, p0=[ 1, 1, 1, 1 ] )
print( popt )
print( pcov ) 
### always look at pcov
### shows that fit is basically insenitive to parameter D

xfull = np.linspace( -3, 15, 155 )
yfull = objective( xfull, *popt )

fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )

ax.plot( xl, yl, ls='', marker='o' )
ax.plot( xfull, yfull )

plt.show()

请注意，通过*导入所有内容通常不是一个好主意。