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

Posted 2023-03-12

【中文标题】使用 Python 和 lmfit 拟合复杂模型？

【英文标题】：Fitting complex model using Python and lmfit?

【发布时间】：2014-06-18 22:00:17

【问题描述】：

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

我可以尝试使用shared parameters 或piecewise 方法将问题分离为实部和虚部，但有没有办法直接使用复杂函数来实现？ 仅拟合函数的实部效果很好，但是当我定义复杂的残差函数时，我得到：

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)
    """
    #=================================================================

编辑： 我解决了虚部和实部分离变量的问题。数据的形状应为[[imaginary_data],[real_data]]，目标函数必须返回一维数组。

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 ...

【问题讨论】：

虽然我不确定，但我认为您得到的 TypeError 基本上回答了您的问题：lmfit 需要知道是增加还是减少其参数，并且由于订购对于未定义的复数，它不能（j 大于还是小于 1？）。

有时答案太明显以至于看不到它...我解决了将探针分离为实变量和虚变量的问题，如果有兴趣，请参阅我的编辑。谢谢埃弗特！ :)

@JurKravla 您能否请answer your own question，而不是在您的问题的编辑中显示解决方案？这个问题不断出现在 SciPy 标签的“未回答的问题”中。

【参考方案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()