lmfit -py使用數組進行參數優化

情況：
我正在嘗試優化天然小溪的參數，其中氣體根據合理建立的方程以特定的速率脫氣或排氣。我們在下游的某個距離處測量濃度，並希望使用優化技術來確定模型中某些未知參數的值。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/

https://github.com/lmfit/lmfit-py

http://cars9.uchicago.edu/software/python/lmfit/

來源

2014-03-05 Oregano

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

，然後使用這些在目標函數中標量「每數據集」參數。當然，您仍然可以對所有數據集使用全局參數。

來源

2014-03-09 13:06:38

lmfit -py使用數組進行參數優化

回答

相關問題