lmfit -py 使用数组进行参数优化
情况：
我正在尝试优化一个自然溪流的参数，在这个溪流中，气体会以一定的速率释放或吸收，这个速率是根据一些比较成熟的公式来计算的。我们在溪流下游的某些距离上测量了气体的浓度，现在想用优化技术来找出模型中一些未知参数的值。
我正在使用 lmfit lmfit-py (github) 来优化参数，下面是我用的代码。
我希望某些未知参数可以在溪流的不同采样点上有不同的值（这更符合实际情况）。在脚本中，这些参数被称为“局部参数”。因此，这些参数的结果应该是每个位置的一个列表。
其他参数应该优化为在溪流的所有位置上都有相同的值，这被称为“全局参数”。因此，这些参数的优化结果应该是一个单一的值。
我使用的是 Python 3.2。
问题：
- 目前我得到的所有参数的结果都是单一值。
lmfit能否用来生成优化后的数组？- 我觉得我可能没有正确设置脚本，可能没有正确调用边界和初始猜测，或者其他什么问题？？？
我使用的脚本：
import numpy as np
import matplotlib.pyplot as plt
from lmfit import minimize, Parameters, report_fit
import scipy as sp
# import data to be fitted
x,dc_dx,E,lmbd,c,h,th,theta,d,w,c_min,\
c_max,Q_min,Q_max,w_min,w_max,d_min,d_max,theta_min,theta_max,\
I_min,I_max,k_min,k_max,h_min,h_max,th_min,th_max,c_ini,Q_ini,\
w_ini,d_ini,theta_ini,I_ini,k_ini,h_ini,th_ini,cigl_min,\
cigl_max,gammagl_min,gammagl_max,\
cigl_ini,gammagl_ini,kgl_min,kgl_max,kgl_ini,\
temp,scond,econd,disol,O2per,O2,pH,ORP,\
c_atmdef,dO2_dx,kO2_ini,kO2_min,kO2_max,O2i=\
np.genfromtxt("Rn2011DryInput_1.csv",unpack=True,\
skip_header =2,delimiter=",")
'''
'Global parameters' are those which when optimised result in the
same value along the whole transect.
'Local parameters' are those parameters optimised where the value can
vary at each sampling site.
(global parameters should be packed at the beginning)
for globals unpack each one e.g. p1=par[0]
p2 = par[1]
'''
# Define smoothing error
def ratio(a,b):
error = 0.0
if b > sys.float_info.epsilon:
error = ((a-b)/b)*2
else:
if a > sys.float_info.epsilon:
error = ((a-b)/a)**2
# Try reducing the error for the regularisation
error = 0.5*error
return(error)
# Define model
def radon_f(par, *args):
nargs = 12 # number of arguments
npl = 4 # number of local parameters
npg = 2 # number of global parametrs
ci=par['ci'].value
gamma=par['gamma'].value
n = int(len(par)/npl)
error = 0.0
for i in range(n):
Q=par['Q'].value
I=par['I'].value
k=par['k'].value
kO2=par['kO2'].value
dc_dx=args[i*nargs]
E=args[i*nargs+1]
lmbd=args[i*nargs+2]
c=args[i*nargs+3]
h=args[i*nargs+4]
th=args[i*nargs+5]
theta=args[i*nargs+6]
d=args[i*nargs+7]
w=args[i*nargs+8]
O2i=args[i*nargs+9]
O2=args[i*nargs+10]
c_atmdef=args[i*nargs+10]
measurements = dc_dx
model=(I*(ci-c)+w*E*c-k*w*c-d*w*lmbd*c+\
gamma*h*w*theta/(1+lmbd*th)-lmbd*h*w*theta*c/(1+lmbd*th))/Q
error= error + ((measurements - model)/measurements)**2
measurements=dO2_dx
model = (I*(O2i-O2)+w*E*O2-kO2*w*c_atmdef)/Q
error= error + ((measurements - model)/measurements)**2
if i > 0: # apply regularization
Qp,Ip,kp,kO2p=par[(npg+npl*(i-1)):(npg+npl*i)]
error = error + ratio(Q,Qp)
error = error + ratio(I,Ip)
error = error + ratio(k,kp)
error = error + ratio(kO2,kO2p)
return(error)
# create a set of parameters
par = Parameters()
for i in range(int(len(x))):
par.add('Q', value=Q_ini[i], vary=True, min= Q_min[i], max=Q_max[i])
par.add('I', value=I_ini[i], vary=True, min= I_min[i], max=I_max[i])
par.add('k', value=k_ini[i], vary=True, min= 0.1, max=6.2 )
par.add('ci', value=cigl_ini[i], vary=False, min= 6000, max=25000 )
par.add('gamma',value=gammagl_ini[i], vary=False, min=200, max=3000 )
par.add('kO2', value=kO2_ini[i], vary=True, min=2.5, max=10 )
# Define arguements
args=[]
for i in range(len(E)):
args.append(dc_dx[i])
args.append(E[i])
args.append(lmbd[i])
args.append(c[i])
args.append(h[i])
args.append(th[i])
args.append(theta[i])
args.append(d[i])
args.append(w[i])
args.append(O2i[i])
# run the fit
result = minimize(radon_f, par, args=args)
print("all results")
print(result)
Q=par['Q'].value
I=par['I'].value
k=par['k'].value
kO2=par['kO2'].value
print(Q)
print(I)
print(k)
print(kO2)
# print results of global parameters
print('ci', 'gamma')
print(result[0][0], result[0][3])
for i in range(len(c_min)):
# Print out the results of local parameters
print(result[0][npg+npl*i],result[0][npg+npl*i+1],\
result[0][npg+npl*i+2],'\t')
我使用的数据：
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
x dc_dx E lmbd c h th theta d w c_min c_max Q_min Q_max w_min w_max d_min d_max theta_min theta_max I_min I_max k_min k_max h_min h_max th_min th_max c_ini Q_ini w_ini d_ini theta_ini I_ini k_ini h_ini th_ini cigl_min cigl_max gammagl_min gammagl_max cigl_ini gammagl_ini kgl_min kgl_max kgl_ini temp scond econd disol O2per O2 pH ORP c_atmdef dO2_dx kO2_ini kO2_min kO2_max O2i
514 -6.763285345 0.0045 0.181 9572 0.7 0.5 0.3 0.5 12.2 8614.762109 10529.15369 0 14200 3 15 0.1 2 0.2 0.45 0 1000 1 6.17 0.4 1 0 2 9572 14000 3 0.4 0.3 50 2 0.6 0.6 6000 25000 200 3000 7000 700 1 7 2 26.93 752.3 780.1 489 43 3.42 7.27 8 -267.48 0.012840237 5 2.5 10 1.3
683 -5.490998291 0.0045 0.181 8429 0.7 0.5 0.3 0.5 12.2 7586.066408 9271.858943 0 14200 3 15 0.1 2 0.2 0.45 0 1000 1 6.17 0.4 1 0 2 8429 14000 3 0.4 0.3 50 2 0.6 0.6 26.61 750.4 773.4 487.7 69.8 5.59 7.33 12 -265.31 0.005837412 5 2.5 10 1.3
949 -4.22793979 0.0045 0.181 7307 0.7 0.5 0.3 0.5 12.2 6576.106938 8037.464035 0 14200 3 15 0.1 2 0.2 0.45 0 1000 1 6.17 0.4 1 0 2 7307 14000 3 0.4 0.3 50 2 0.6 0.6 26.3 750.7 769.4 488 65.6 5.28 7.38 26 -265.62 0.000472201 5 2.5 10 1.3
1287 -2.085993575 0.0045 0.181 5875 0.7 0.5 0.3 0.5 12.2 5287.160328 6462.084846 0 14200 3 15 0.1 2 0.2 0.45 0 1000 1 6.17 0.4 1 0 2 5875 14000 3 0.4 0.3 50 2 0.6 0.6 26.15 750 766 487 71.6 5.78 7.47 19.5 -265.12 0.00183869 5 2.5 10 1.3
1623 -1.048348565 0.0045 0.181 5897 0.7 0.5 0.3 0.5 12.2 5306.871121 6486.175814 12822.3 15671.7 3 15 0.1 2 0.2 0.45 0 1000 1 6.17 0.4 1 0 2 5897 14247 3 0.4 0.3 50 2 0.6 0.6 25.98 748.1 762.2 486.3 80.5 6.52 7.64 27 -264.38 0.001575445 5 2.5 10 1.3
1992 3.150847718 0.0045 0.181 5099 0.7 0.5 0.3 0.5 12.2 4588.91133 5608.669404 14500 16500 3 15 0.1 2 0.2 0.45 0 1000 1 6.17 0.4 1 0 2 5099 16707 8 1 0.3 50 2 0.6 0.6 26.03 750 765 487 85 6.87 7.56 26 -264.03 -0.003467278 5 2.5 10 1.3
2488 17.92239204 0.0045 0.181 9297 0.7 0.5 0.3 0.5 12.2 8367.050656 10226.39525 35500 59500 3 15 0.1 2 0.2 0.45 0 1000 1 6.17 0.4 1 0 2 9297 49279 8 1 0.3 100 2 0.6 0.6 29.06 726 783 472 38.5 2.96 7.14 -83.4 -267.94 -0.010477451 5 2.5 10 1.3
2569 10.03067057 0.0045 0.181 11515 0.7 0.5 0.3 0.5 12.2 10363.14089 12666.06109 44500 60000 3 15 0.1 2 0.2 0.45 0 1000 1 6.17 0.4 1 0 2 11515 59739 8 1 0.3 750 2 0.6 0.6 29.5 730 793 474 43 3.28 7.07 -1.3 -267.62 0.002783579 5 2.5 10 1.3
2835 -6.606394762 0.0045 0.181 9568 0.7 0.5 0.3 0.5 12.2 8610.764203 10524.26736 44500 60000 3 15 0.1 2 0.2 0.45 0 1000 1 6.17 0.4 1 0 2 9568 58576 8 1 0.3 250 2 0.6 0.6 29.3 727 787 472 48 3.71 7.12 11.1 -267.19 0.001373823 5 2.5 10 1.3
3224 -4.694876493 0.0045 0.181 7275 0.7 0.5 0.3 0.5 12.2 6547.652797 8002.686751 44500 60000 3 15 0.1 2 0.2 0.45 0 1000 1 6.17 0.4 1 0 2 7275 56874 8 1 0.3 150 2 0.6 0.6 28.29 729.9 775.7 474.4 53.4 4.15 7.27 21 -266.75 0.001830971 5 2.5 10 1.3
3609 -4.654680461 0.0045 0.181 5929 0.7 0.5 0.3 0.5 12.2 5336.000281 6521.778121 44500 60000 3 15 0.1 2 0.2 0.45 0 1000 1 6.17 0.4 1 0 2 5929 55190 10 1 0.3 150 2 0.6 0.6 29.403142 734 702 477 59.6 5.13 7.22 -28 -265.77 0.001991543 5 2.5 10 1.3
4082 -3.263752353 0.0045 0.181 3180 0.7 0.5 0.3 0.5 12.2 2861.606999 3497.519665 44500 60000 3 15 0.1 2 0.2 0.45 0 1000 1 6.17 0.4 1 0 2 3180 53121 10 1.5 0.3 150 2 0.6 0.6 28.164949 729 681 474 66 5.81 7.5 -33 -265.09 0.001000506 5 2.5 10 1.3
5218 1.167013246 0.0045 0.181 2367 0.7 0.5 0.3 0.5 12.2 2130.615084 2604.085102 44500 60000 3 15 0.1 2 0.2 0.45 0 1000 1 6.17 0.4 1 0 2 2367 48152 10 1.5 0.3 50 2 0.6 0.6 26.425339 684 616 444 70.8 6.45 7.83 -41 -264.45 0.00048454 5 2.5 10 1.3
5506 1.045825103 0.0045 0.181 3245 0.7 0.5 0.3 0.5 12.2 2920.916644 3570.009232 44500 60000 3 15 0.1 2 0.2 0.45 0 1000 0 6.17 0.4 1 0 2 3245 46893 10 1.5 0.3 150 2 0.6 0.6 26.527669 681 615 443 75.4 6.85 7.83 -52 -264.05 0.000191651 5 2.5 10 1.3
6043 -0.741605916 0.0045 0.181 2731 0.7 0.5 0.3 0.5 12.2 2458.228071 3004.500975 40089.6 48998.4 3 15 0.1 2 0.2 0.45 0 1000 0 6.17 0.4 1 0 2 2731 44544 10 1.5 0.3 100 2 0.6 0.6 24.900622 690 602 448 67.2 6.31 7.75 -38 -264.59 0.000307351 5 2.5 10 1.3
7851 -0.366090159 0.0045 0.181 1781 0.7 0.5 0.3 0.5 12.2 1602.550136 1958.672388 44500 53500 3 15 0.1 2 0.2 0.45 0 1000 0 6.17 0.4 1 0 2 1781 46298 10 1.5 0.3 50 2 0.6 0.6 24.675496 679 590 441 71.8 6.77 7.85 -26 -264.13 0.000430647 5 2.5 10 1.3
15277 -0.206321214 0.0045 0.181 248 0.7 0.5 0.3 0.5 12.2 223.6229375 273.3169236 48150 58850 3 15 0.1 2 0.2 0.45 0 1000 0 6.17 0.4 1 0 2 248 53500 10 1.5 0.3 25 2 0.6 0.6 22.97 621.9 597.8 404.2 114.9 9.84 8.02 11 -261.06 0.000413412 5 2.5 10 1.3
非常感谢！
相关背景信息：
http://lmfit.github.io/lmfit-py/
https://pypi.python.org/pypi/lmfit/
1 个回答
1
LMFIT 这个工具不能直接优化一组值。如果你想对多个“数据集”进行拟合，唯一的方法就是自己写一个目标函数来实现这一点——看起来你已经有这个了。不过，你在创建参数时用的
# create a set of parameters
par = Parameters()
for i in range(int(len(x))):
par.add('Q', value=Q_ini[i], vary=True, min= Q_min[i], max=Q_max[i])
par.add('I', value=I_ini[i], vary=True, min= I_min[i], max=I_max[i])
....
其实是在覆盖参数 Q 和 I 的定义。你可能想要为每个数据集单独设置这些参数，可以用类似下面的方式
# create a set of parameters
par = Parameters()
for i in range(int(len(x))):
par.add('Q_%i'%(i+1), value=Q_ini[i], vary=True, min= Q_min[i], max=Q_max[i])
par.add('I_%i'%(i+1), value=I_ini[i], vary=True, min= I_min[i], max=I_max[i])
....
然后在你的目标函数中使用这些针对每个数据集的参数。当然，你也可以有一些参数是对所有数据集都适用的。