sir模型python拟合_使用Python和lmfit拟合复杂的模型？

I would like to fit ellipsometric data to complex model using LMFit. Two measured parameters, psi and delta, are variables in a complex function rho.

I could try with separating problem to real and imaginary part with shared parameters or piecewise approach, but is there any way to do it directly with complex function?

Fitting only real part of function works beautifully, but when I define complex residual function I get:

TypeError: no ordering relation is defined for complex numbers.

Below is my code for real function fitting and my attempt at tackling complex fit problem:

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)

"""

#=================================================================

Edit:

I solved problem with separated variables for imaginary and real part. Data should be shaped as [[imaginary_data],[real_data]], objective function must return 1D array.

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

解决方案

The lmfit FAQ suggests simply taking both real and imaginary parts by using numpy.ndarray.view, which means you don't need to go through the separation of the real and imaginary parts manually.

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()