6

我想使用 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)
    """
    #=================================================================

编辑： 我解决了虚部和实部分离变量的问题。数据的形状应为 [[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 ...
4

1 回答 1

1

lmfit 常见问题解答建议通过using 简单地同时获取实部和虚部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()
于 2018-06-09T08:23:11.677 回答