状況:
十分に確立された方程式に従って、ガスが一定の速度で脱ガスまたはガスを吸入する自然の小川のパラメータを最適化しようとしています。下流の特定の距離で測定された濃度があり、最適化手法を使用して、モデル内のいくつかの未知のパラメーターの値を決定したいと考えています。
以下に貼り付けたコードを使用して、 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/