使用Python和lmfit拟合复杂模型？

Question

我想使用LMFit将椭圆测量数据拟合到复杂模型中。两个测量参数psi和delta是复函数rho的变量。

我可以尝试使用共享参数或分段方法，将问题分解为实部和虚部。但是否有直接使用复函数进行处理的方法？

仅拟合函数的实部效果良好，但当我定义复残差函数时，会出现以下错误：

TypeError：未为复数定义排序关系。

以下是我用于实函数拟合的代码以及尝试解决复杂拟合问题的代码：

    from __future__ import division
    from __future__ import print_function
    import numpy as np
    from pylab import *
    from lmfit import minimize, Parameters, Parameter, report_errors


    #=================================================================
    #             MODEL

    def r01_p(eps2, th):
        c=cos(th)
        s=(sin(th))**2

        stev= sqrt(eps2) * c - sqrt(1-(s / eps2))
        imen= sqrt(eps2) * c + sqrt(1-(s / eps2))
        return stev/imen

    def r01_s(eps2, th):
        c=cos(th)
        s=(sin(th))**2

        stev= c - sqrt(eps2) * sqrt(1-(s/eps2))
        imen= c + sqrt(eps2) * sqrt(1-(s/eps2))
        return stev/imen


    def rho(eps2, th):
        return r01_p(eps2, th)/r01_s(eps2, th)

    def psi(eps2, th):
        x1=abs(r01_p(eps2, th))
        x2=abs(r01_s(eps2, th))
        return np.arctan2(x1,x2)

    #=================================================================
    #                   REAL FIT
    #

    #%%

    # generate data from model  
    th=linspace(deg2rad(45),deg2rad(70),70-45)
    error=0.01
    var_re=np.random.normal(size=len(th), scale=error)
    data = psi(2,th) + var_re

    # residual function
    def residuals(params, th, data):
        eps2 = params['eps2'].value

        diff = psi(eps2, th) - data
        return diff

    # create a set of Parameters
    params = Parameters()
    params.add('eps2',   value= 1.0,  min=1.5, max=3.0)


    # do fit, here with leastsq model
    result = minimize(residuals, params, args=(th, data),method="leastsq")

    # calculate final result
    final = data + result.residual

    # write error report
    report_errors(params)


    # try to plot results
    th, data, final=rad2deg([th, data, final])
    try:
        import pylab
        clf()
        fig=plot(th, data, 'r o',
                 th, final, 'b')
        setp(fig,lw=2.)
        xlabel(r'$\theta$ $(^{\circ})$', size=20)
        ylabel(r'$\psi$ $(^{\circ})$',size=20)

        show()

    except:
        pass

    #%%
    #=================================================================
    #                   COMPLEX FIT

    # TypeError: no ordering relation is defined for complex numbers

    """
    # data from model with added noise   
    th=linspace(deg2rad(45),deg2rad(70),70-45)
    error=0.001
    var_re=np.random.normal(size=len(th), scale=error)
    var_im=np.random.normal(size=len(th), scale=error) * 1j

    data = rho(4-1j,th) + var_re + var_im


    # residual function
    def residuals(params, th, data):
        eps2 = params['eps2'].value

        diff = rho(eps2, th) - data
        return np.abs(diff)

    # create a set of Parameters
    params = Parameters()
    params.add('eps2',   value= 1.5+1j,  min=1+1j, max=3+3j)


    # do fit, here with leastsq model
    result = minimize(residuals, params, args=(th, data),method="leastsq")

    # calculate final result
    final = data + result.residual

    # write error report
    report_errors(params)
    """
    #=================================================================

编辑： 我通过将虚部和实部分开解决了问题。数据应该被形状化为[[虚数数据],[实数数据]]，目标函数必须返回一维数组。

def objective(params, th_data, data):
    eps_re  = params['eps_re'].value
    eps_im  = params['eps_im'].value
    d       = params['d'].value

    residual_delta = data[0,:] - delta(eps_re - eps_im*1j, d, frac, lambd, th_data)
    residual_psi   = data[1,:] - psi(eps_re - eps_im*1j, d, frac, lambd, th_data)

    return np.append(residual_delta,residual_psi)

# create a set of Parameters
params = Parameters()
params.add('eps_re',   value= 1.5,  min=1.0,      max=5  )
params.add('eps_im',   value= 1.0,  min=0.0,      max=5  )
params.add('d',        value= 10.0,   min=5.0,    max=100.0   )


# All available methods
methods=['leastsq','nelder','lbfgsb','anneal','powell','cobyla','slsqp']
# Chosen method
#metoda='leastsq'

# run the global fit to all the data sets
result = minimize(objective, params, args=(th_data,data),method=metoda))

....

return ...

- Jur Kravla

---

虽然我不确定，但我认为你得到的“TypeError”基本上回答了你的问题：lmfit需要知道是增加还是减少其参数，由于复数的排序未定义，因此它无法进行操作（j大于还是小于1？）。 - user707650

有时候答案太明显了，以至于我们看不到它... 我通过将问题分离为实数和虚数变量来解决问题，如果感兴趣可以查看我的编辑。 谢谢Evert！ :) - Jur Kravla

@JurKravla，您能否请自行回答您自己的问题，而不是在问题的编辑中展示解决方案？这个问题一直出现在SciPy标签的“未回答的问题”中。 - chthonicdaemon

Answer 1

lmfit FAQ建议使用numpy.ndarray.view来简单地取实部和虚部，这意味着您不需要手动进行实部和虚部的分离。

def residuals(params, th, data):
    eps2 = params['eps2'].value

    diff = rho(eps2, th) - data
    # The only change required is to use view instead of abs.
    return diff.view()

- chthonicdaemon