问题描述
我想使用LMFit使椭偏数据适合复杂模型.测得的两个参数psi
和delta
是复杂函数rho
中的变量.
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:
下面是我的实函数拟合代码,也是我尝试解决复杂拟合问题的代码:
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)
"""
#=================================================================
修改:我用虚部和实部分开的变量解决了问题.数据的形状应为[[imaginary_data],[real_data]],目标函数必须返回一维数组.
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 ...
推荐答案
lmfit常见问题解答建议使用 numpy.ndarray.view
,这意味着您无需手动进行实部和虚部的分离.
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()
这篇关于使用Python和lmfit拟合复杂的模型?的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!