0

情况：
我正在尝试根据合理建立的方程式优化天然小溪的参数，其中气体以一定的速率脱气或加气。我们在下游特定距离处测量了浓度，并希望使用优化技术来确定模型中某些未知参数的值。

我正在尝试使用 lmfit lmfit-py (github)使用下面粘贴的代码优化参数。

  • 我希望允许一些未知参数对于流中的每个不同采样点具有不同的值（这反映了实际情况）。在脚本中，这些被称为“本地参数”。因此，这些参数的结果应该是流中每个位置的列表。

  • 其他参数应优化为沿流的所有位置具有相同的值，称为“全局参数”。因此，这些参数的优化结果应该是单个值。

我正在使用 Python 3.2。

问题：

  1. 目前我只得到所有参数的单值结果。
  2. 可以lmfit用来产生优化的数组吗？
  3. 我真的认为我没有正确设置脚本，可能没有正确调用边界和初始猜测或其他什么？？？

我正在使用的脚本： 

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/

4

1 回答 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])
    ....

然后在目标函数中使用那些标量“每个数据集”参数。当然，您仍然可以为所有数据集设置全局参数。

于 2014-03-09T13:06:38.570 回答