我正在使用 lmfit 来寻找优化模型与分子光谱的拟合的参数。残差程序调用一个 Fortran 代码,该代码计算能级并将它们与可用的实验数据匹配,提供作为 lmfit-minimize 输入的残差集。
有两种可能的情况:线性分子和非线性分子。我为他们每个人准备了一个例子。使用线性分子,输出是我所期望的,达到与我以前使用 Minuit 的 Fortran 版本进行参数优化的代码相同的最小值。这是一个示例,其中打印的值是每次调用函数时获得的卡方值。
import numpy as np
from lmfit import minimize, Parameters, fit_report
from subprocess import Popen, PIPE
#
%run ../../bin/residuals_U3.py
input_filename = "HNC_input_file"
BENT = ".F."
exp_data_file="exp_HNC_Danielle.dat"
N_val=60
LMAX=7
VMAX=7
EMINL=".F."
U3H_FIT_param=Parameters()
U3H_FIT_param.add('P11',value=2390.0, vary=True)
U3H_FIT_param.add('P21',value=-37.0,vary=True)
U3H_FIT_param.add('P22',value=21.0, vary=True)
U3H_FIT_param.add('P23',value=-9.0, vary=True)
out_0 = minimize(residuals_U3, U3H_FIT_param, args=(input_filename, BENT,
exp_data_file, N_val, LMAX, VMAX, EMINL))
31298390.0981
31298390.0981
31298390.0981
31298405.9262
31298393.7425
31298390.2979
31298399.7958
1239848240.71
....
5.7028931985
5.70223335575
5.68090132895
5.67973533085
5.68050720657
5.68089119394
5.68021995958
5.57489186835
5.57489701562
5.57489288067
5.57489188595
5.57489412245
5.57489158478
print fit_report(out_0)
[[Fit Statistics]]
# function evals = 169
# data points = 19
# variables = 4
chi-square = 5.575
reduced chi-square = 0.372
Akaike info crit = -15.297
Bayesian info crit = -11.519
[[Variables]]
P11: 2375.84727 +/- 59.96417 (2.52%) (init= 2390)
P21: -37.5581932 +/- 1.065396 (2.84%) (init=-37)
P22: 19.2964371 +/- 0.426936 (2.21%) (init= 21)
P23: -9.63266677 +/- 0.261630 (2.72%) (init=-9)
P31: 0 (fixed)
P32: 0 (fixed)
P33: 0 (fixed)
P41: 0 (fixed)
P42: 0 (fixed)
P43: 0 (fixed)
P44: 0 (fixed)
P45: 0 (fixed)
P46: 0 (fixed)
P47: 0 (fixed)
[[Correlations]] (unreported correlations are < 0.100)
C(P11, P23) = -1.000
C(P11, P21) = -1.000
C(P21, P23) = 1.000
C(P22, P23) = -1.000
C(P11, P22) = 1.000
C(P21, P22) = -1.000 code here
到目前为止,一切都很好。一旦我尝试拟合来自非线性分子种类的数据,就会出现问题。然后,经过很少的评估,程序似乎收敛到提供的初始值。如果我们将此示例称为 1:
import numpy as np
from lmfit import minimize, Parameters, fit_report
from subprocess import Popen, PIPE
%run ../../bin/residuals_U3.py
input_filename = "NCNCS_input_file"
BENT = ".T."
exp_data_file="test_exp_NCNCS.dat"
N_val=70
LMAX=9
VMAX=6
EMINL=".F."
# Parameters far from the minimum I already know
U3H_FIT_param=Parameters()
U3H_FIT_param.add('P11',value=200.0, vary=True)
U3H_FIT_param.add('P21',value=-2.60,vary=True)
U3H_FIT_param.add('P22',value=1.50, vary=True)
U3H_FIT_param.add('P23',value=-0.8, vary=True)
out_0 = minimize(residuals_U3, U3H_FIT_param, args=(input_filename, BENT, exp_data_file, N_val, LMAX, VMAX, EMINL))
33820.0608635
33820.0608635
33820.0608635
33942.2383744
34935.7037896
33810.3126988
33357.9527561
33303.3818229
34574.0192135
34526.2057525
33337.2494285
34458.8982496
34798.9915284
print fit_report(out_0)
[[Fit Statistics]]
# function evals = 13
# data points = 70
# variables = 4
chi-square = 33303.382
reduced chi-square = 504.597
Akaike info crit = 439.544
Bayesian info crit = 448.538
[[Variables]]
P11: 199.999988 +/- 4.02e-05 (0.00%) (init= 200)
P21: -2.60000161 +/- 6.71e-07 (0.00%) (init=-2.6)
P22: 1.50000025 +/- 3.57e-07 (0.00%) (init= 1.5)
P23: -0.79999987 +/- 2.12e-07 (0.00%) (init=-0.8)
P31: 0 (fixed)
P32: 0 (fixed)
P33: 0 (fixed)
P41: 0 (fixed)
P42: 0 (fixed)
P43: 0 (fixed)
P44: 0 (fixed)
P45: 0 (fixed)
P46: 0 (fixed)
P47: 0 (fixed)
[[Correlations]] (unreported correlations are < 0.100)
C(P21, P23) = -0.675
C(P21, P22) = 0.506
C(P22, P23) = -0.430
C(P11, P21) = -0.420
C(P11, P23) = -0.317
C(P11, P22) = -0.216
如果我们使用更详细的输出重复示例 1,除了之前的 Chi^2 之外,还包括我们获得的参数值和残差:
U3H_FIT_param=Parameters()
U3H_FIT_param.add('P11',value=200.0, vary=True)
U3H_FIT_param.add('P21',value=-2.60,vary=True)
U3H_FIT_param.add('P22',value=1.50, vary=True)
U3H_FIT_param.add('P23',value=-0.8, vary=True)
out_0 = minimize(residuals_U3, U3H_FIT_param, args=(input_filename, BENT, exp_data_file, N_val, LMAX, VMAX, EMINL))
&INP1b P11=200.0 /
&INP2b P21=-2.6, P22=1.5, P23=-0.8 /
RESIDUALS [ 0. -3.75658039 -10.85664646 -18.29270568 -21.23633151
-23.97010961 -27.99053887 -0.28486237 -3.86773685 -9.31351971
-15.92856641 -20.2686777 -24.7703342 -30.12984172 -0.8655828
-4.27628737 -9.23390995 -15.22920459 -20.20446745 -25.50896835
-31.54014503 -1.98904718 -5.28870374 -9.64523663 -14.86388301
-20.30393994 -26.28976779 -32.87036748 -3.4154388 -5.99093571
-10.52895266 -15.51556868 -21.64134653 -27.44348466 -35.0922282
-4.59343935 -7.18069815 -11.05888025 -16.24498349 -22.25796287
-29.17718901 -37.86040995 -6.82393206 -8.37059213 -12.37991401
-17.32380238 -23.6259757 -30.77116529 -39.6603254 -7.37663293
-9.4855129 -13.15066135 -18.13207261 -24.54463108 -33.16897104
-43.02725922 -8.32876467 -10.8207887 -14.54627353 -19.60041984
-26.69368656 -34.95192225 -46.08130388 -9.44803578 -11.71293473
-15.45927873 -21.18174524 -28.33172007 -37.40255102 -49.41886369]
33820.0608635
&INP1b P11=200.0 /
&INP2b P21=-2.6, P22=1.5, P23=-0.8 /
RESIDUALS [ 0. -3.75658039 -10.85664646 -18.29270568 -21.23633151
-23.97010961 -27.99053887 -0.28486237 -3.86773685 -9.31351971
-15.92856641 -20.2686777 -24.7703342 -30.12984172 -0.8655828
-4.27628737 -9.23390995 -15.22920459 -20.20446745 -25.50896835
-31.54014503 -1.98904718 -5.28870374 -9.64523663 -14.86388301
-20.30393994 -26.28976779 -32.87036748 -3.4154388 -5.99093571
-10.52895266 -15.51556868 -21.64134653 -27.44348466 -35.0922282
-4.59343935 -7.18069815 -11.05888025 -16.24498349 -22.25796287
-29.17718901 -37.86040995 -6.82393206 -8.37059213 -12.37991401
-17.32380238 -23.6259757 -30.77116529 -39.6603254 -7.37663293
-9.4855129 -13.15066135 -18.13207261 -24.54463108 -33.16897104
-43.02725922 -8.32876467 -10.8207887 -14.54627353 -19.60041984
-26.69368656 -34.95192225 -46.08130388 -9.44803578 -11.71293473
-15.45927873 -21.18174524 -28.33172007 -37.40255102 -49.41886369]
33820.0608635
&INP1b P11=200.0 /
&INP2b P21=-2.6, P22=1.5, P23=-0.8 /
RESIDUALS [ 0. -3.75658039 -10.85664646 -18.29270568 -21.23633151
-23.97010961 -27.99053887 -0.28486237 -3.86773685 -9.31351971
-15.92856641 -20.2686777 -24.7703342 -30.12984172 -0.8655828
-4.27628737 -9.23390995 -15.22920459 -20.20446745 -25.50896835
-31.54014503 -1.98904718 -5.28870374 -9.64523663 -14.86388301
-20.30393994 -26.28976779 -32.87036748 -3.4154388 -5.99093571
-10.52895266 -15.51556868 -21.64134653 -27.44348466 -35.0922282
-4.59343935 -7.18069815 -11.05888025 -16.24498349 -22.25796287
-29.17718901 -37.86040995 -6.82393206 -8.37059213 -12.37991401
-17.32380238 -23.6259757 -30.77116529 -39.6603254 -7.37663293
-9.4855129 -13.15066135 -18.13207261 -24.54463108 -33.16897104
-43.02725922 -8.32876467 -10.8207887 -14.54627353 -19.60041984
-26.69368656 -34.95192225 -46.08130388 -9.44803578 -11.71293473
-15.45927873 -21.18174524 -28.33172007 -37.40255102 -49.41886369]
33820.0608635
.....
&INP1b P11=199.999988473 /
&INP2b P21=-2.60000161305, P22=1.50000025484, P23=-0.799999875067 /
RESIDUALS [ 0. -3.99268261 -11.73264434 -19.05089687 -21.8743998
-23.98786565 -28.05424445 -0.80354366 -4.29117203 -10.41169177
-16.11475578 -20.76739082 -24.72907688 -29.96025254 -1.19554552
-4.77165975 -10.04051152 -15.40288703 -20.18713017 -25.8307694
-31.91088782 -2.46101195 -5.74047469 -9.9850673 -15.24025329
-20.72837057 -26.68827531 -33.31621124 -3.66948178 -6.53021075
-10.89700164 -15.81327816 -21.81723225 -28.08793728 -35.76616676
-4.99405722 -7.65935238 -11.55129453 -16.5529505 -22.40531267
-29.36326533 -38.30271257 -7.21713077 -8.76345313 -12.77142684
-17.71182548 -24.01827081 -31.19676872 -40.10579542 -7.61666334
-9.92207061 -13.58164103 -18.52451949 -25.03610104 -33.39256778
-43.32754575 -8.67692476 -11.12324272 -15.00726898 -20.04660062
-27.02514438 -35.41059799 -46.51995156 -10.18456526 -12.0963503
-15.60152462 -21.17712212 -28.94048104 -37.94146875 -49.83299737]
34798.9915284
print fit_report(out_0)
[[Fit Statistics]]
# function evals = 13
# data points = 70
# variables = 4
chi-square = 33303.382
reduced chi-square = 504.597
Akaike info crit = 439.544
Bayesian info crit = 448.538
[[Variables]]
P11: 199.999988 +/- 4.02e-05 (0.00%) (init= 200)
P21: -2.60000161 +/- 6.71e-07 (0.00%) (init=-2.6)
P22: 1.50000025 +/- 3.57e-07 (0.00%) (init= 1.5)
P23: -0.79999987 +/- 2.12e-07 (0.00%) (init=-0.8)
P31: 0 (fixed)
P32: 0 (fixed)
P33: 0 (fixed)
P41: 0 (fixed)
P42: 0 (fixed)
P43: 0 (fixed)
P44: 0 (fixed)
P45: 0 (fixed)
P46: 0 (fixed)
P47: 0 (fixed)
[[Correlations]] (unreported correlations are < 0.100)
C(P21, P23) = -0.675
C(P21, P22) = 0.506
C(P22, P23) = -0.430
C(P11, P21) = -0.420
C(P11, P23) = -0.317
C(P11, P22) = -0.216
如示例(示例 2)所示,显然有一个更好的最小值:
# Parameters closer to the minimum
U3H_FIT_param=Parameters()
U3H_FIT_param.add('P11',value=201.0, vary=True)
U3H_FIT_param.add('P21',value=-2.570,vary=True)
U3H_FIT_param.add('P22',value=1.4980, vary=True)
U3H_FIT_param.add('P23',value=-0.81, vary=True)
out_0 = minimize(residuals_U3, U3H_FIT_param, args=(input_filename, BENT, exp_data_file, N_val, LMAX, VMAX, EMINL))
866.43578979
866.43578979
866.43578979
880.37396166
890.227624754
852.194331523
780.074380341
801.333830227
print fit_report(out_0)
[[Fit Statistics]]
# function evals = 8
# data points = 70
# variables = 4
chi-square = 801.334
reduced chi-square = 12.141
Akaike info crit = 178.645
Bayesian info crit = 187.639
[[Variables]]
P11: 200.999990 +/- 7.47e-06 (0.00%) (init= 201)
P21: -2.57000013 +/- 1.08e-07 (0.00%) (init=-2.57)
P22: 1.49800011 +/- 5.88e-08 (0.00%) (init= 1.498)
P23: -0.80999997 +/- 8.11e-09 (0.00%) (init=-0.81)
P31: 0 (fixed)
P32: 0 (fixed)
P33: 0 (fixed)
P41: 0 (fixed)
P42: 0 (fixed)
P43: 0 (fixed)
P44: 0 (fixed)
P45: 0 (fixed)
P46: 0 (fixed)
P47: 0 (fixed)
[[Correlations]] (unreported correlations are < 0.100)
C(P11, P21) = -0.386
C(P11, P22) = -0.379
C(P21, P23) = 0.263
C(P11, P23) = -0.125
但是,如果在这种弯曲的情况下,我将最小化方法更改为 nelder,则结果与预期的一样,即使我开始最小化距离示例 1 中的最小值很远:
U3H_FIT_param=Parameters()
U3H_FIT_param.add('P11',value=200.0, vary=True)
U3H_FIT_param.add('P21',value=-2.60,vary=True)
U3H_FIT_param.add('P22',value=1.50, vary=True)
U3H_FIT_param.add('P23',value=-0.8, vary=True)
out_0 = minimize(residuals_U3, U3H_FIT_param, args=(input_filename, BENT, exp_data_file, N_val, LMAX, VMAX, EMINL), method = "nelder")
33820.0608635
386363.244872
540126.50767
27839.1627683
21180.6786213
398690.128256
105124.802101
98763.3226762
332146.050613
3785.18423466
175555.611803
45570.4112121
15015.0822576
20437.4181311
76313.1360929
3192.21680729
...
355.980240836
332.571174789
332.571174789
355.980240836
350.020618839
358.718935527
332.571174789
332.571174789
332.571174789
332.571174789
print fit_report(out_0)
[[Fit Statistics]]
# function evals = 280
# data points = 70
# variables = 4
chi-square = 332.571
reduced chi-square = 5.039
Akaike info crit = 117.085
Bayesian info crit = 126.079
[[Variables]]
P11: 201.428112 (init= 200)
P21: -2.58733894 (init=-2.6)
P22: 1.52918377 (init= 1.5)
P23: -0.81586079 (init=-0.8)
P31: 0 (fixed)
P32: 0 (fixed)
P33: 0 (fixed)
P41: 0 (fixed)
P42: 0 (fixed)
P43: 0 (fixed)
P44: 0 (fixed)
P45: 0 (fixed)
P46: 0 (fixed)
P47: 0 (fixed)
[[Correlations]] (unreported correlations are < 0.100)
在这种情况下,一些迭代的更详细的输出如下
U3H_FIT_param=Parameters()
U3H_FIT_param.add('P11',value=200.0, vary=True)
U3H_FIT_param.add('P21',value=-2.60,vary=True)
U3H_FIT_param.add('P22',value=1.50, vary=True)
U3H_FIT_param.add('P23',value=-0.8, vary=True)
out_0 = minimize(residuals_U3, U3H_FIT_param, args=(input_filename, BENT, exp_data_file, N_val, LMAX, VMAX, EMINL), method = "nelder")
&INP1b P11=200.0 /
&INP2b P21=-2.6, P22=1.5, P23=-0.8 /
RESIDUALS [ 0. -3.75658039 -10.85664646 -18.29270568 -21.23633151
-23.97010961 -27.99053887 -0.28486237 -3.86773685 -9.31351971
-15.92856641 -20.2686777 -24.7703342 -30.12984172 -0.8655828
-4.27628737 -9.23390995 -15.22920459 -20.20446745 -25.50896835
-31.54014503 -1.98904718 -5.28870374 -9.64523663 -14.86388301
-20.30393994 -26.28976779 -32.87036748 -3.4154388 -5.99093571
-10.52895266 -15.51556868 -21.64134653 -27.44348466 -35.0922282
-4.59343935 -7.18069815 -11.05888025 -16.24498349 -22.25796287
-29.17718901 -37.86040995 -6.82393206 -8.37059213 -12.37991401
-17.32380238 -23.6259757 -30.77116529 -39.6603254 -7.37663293
-9.4855129 -13.15066135 -18.13207261 -24.54463108 -33.16897104
-43.02725922 -8.32876467 -10.8207887 -14.54627353 -19.60041984
-26.69368656 -34.95192225 -46.08130388 -9.44803578 -11.71293473
-15.45927873 -21.18174524 -28.33172007 -37.40255102 -49.41886369]
33820.0608635
&INP1b P11=210.0 /
&INP2b P21=-2.6, P22=1.5, P23=-0.8 /
RESIDUALS [ 0. -26.16300891 -35.92631045 -18.8046449 7.84306164
29.85328814 49.35183801 5.44122598 -7.20663588 -7.51532458
5.28787982 22.79715815 41.85783891 59.23158371 14.12838386
9.2311941 12.82257834 24.05402553 39.05684174 54.70325573
69.54410104 24.02558592 24.14278484 30.06592218 40.05349953
53.44068403 67.14018094 80.79972754 33.95012425 37.30408452
44.53432899 54.88610982 67.23500121 79.31052971 91.39583861
44.20244335 49.60557888 58.08586594 68.72834808 79.93361024
91.24600502 101.35283675 53.69459696 61.88966887 70.56905834
80.97053725 91.89016079 102.16450603 111.61917172 64.64821697
73.47362861 82.86652184 93.10114386 103.25734084 112.65601867
120.69284205 75.35070819 84.29816695 94.28958114 104.26041597
113.8639484 122.29150128 129.26978921 85.42971752 95.71850413
105.5315729 115.31685183 124.24099756 132.05550135 137.78892407]
386363.244872
&INP1b P11=200.0 /
&INP2b P21=-2.73, P22=1.5, P23=-0.8 /
RESIDUALS [ 0. -6.51482403 -11.84722003 -17.94556936 -34.33607138
-57.67632401 -81.40734149 -1.63845707 -8.55295725 -17.17160577
-29.53548333 -47.97384287 -69.24130491 -94.24014416 -4.42106668
-13.37478867 -25.1140549 -40.25940819 -59.31786691 -81.45795179
-106.40406773 -8.47241409 -19.88180645 -33.26697111 -49.98215925
-70.28802414 -93.70400896 -119.18639263 -14.91219834 -26.42593577
-41.64563448 -60.30755635 -81.61793552 -105.65083382 -132.51508599
-20.93474386 -34.42278364 -50.99246305 -70.22964079 -92.3174835
-117.92944692 -146.19607146 -29.37468772 -42.38500645 -59.95420309
-80.42758498 -103.61555313 -130.08637485 -159.51597624 -35.8953851
-51.09068006 -69.33552966 -90.56701074 -114.87032877 -142.75145306
-173.6936809 -43.62317666 -59.85538811 -79.26959621 -101.53354992
-126.78686941 -155.31065489 -187.76051864 -52.18050662 -68.57148933
-88.53569045 -111.87500818 -138.66419476 -168.68410045 -203.08199158]
540126.50767
………………………………………………………………………………
&INP1b P11=201.428112052 /
&INP2b P21=-2.58733894638, P22=1.52918377018, P23=-0.815860796044 /
RESIDUALS [ 0. 2.70481025 2.5647914 -0.66779985 -1.84848475 -2.294670
-2.8511186 -0.06763408 2.30399185 2.4069086 0.60429218 -1.17025016
-1.7678735 -2.25363581 -0.97533799 1.89312414 2.19413958 1.4505511
-0.1054865 -1.2606611 -2.4558844 -1.59327497 0.99634674 1.79137313
1.42330317 0.81319402 -0.59261343 -1.9298181 -2.60320618 0.49113672
1.90950601 1.62687398 1.43875694 0.03109259 -2.37265779 -3.45057765
0.19965851 1.72320248 2.2898178 1.34013775 0.40188506 -2.25368524
-4.79999944 -0.29638363 1.49036912 2.26290378 1.88161992 0.70066741
-2.00055299 -4.2727996 -0.37954819 2.06048067 3.136022 2.91270888
0.68341221 -2.54241374 -4.11982953 -0.25818249 2.16914822 3.46019775
2.91502288 1.04311452 -2.62395441 -3.82618485 0.00951645 2.95514821
3.92300846 3.37200094 1.280833 -3.34939407]
332.571174789
&INP1b P11=201.428112052 /
&INP2b P21=-2.58733894638, P22=1.52918377018, P23=-0.815860796044 /
RESIDUALS [ 0. 2.70481025 2.5647914 -0.66779985 -1.84848475 -2.29467012
-2.8511186 -0.06763408 2.30399185 2.4069086 0.60429218 -1.17025016
-1.7678735 -2.25363581 -0.97533799 1.89312414 2.19413958 1.4505511
-0.1054865 -1.2606611 -2.4558844 -1.59327497 0.99634674 1.79137313
1.42330317 0.81319402 -0.59261343 -1.9298181 -2.60320618 0.49113672
1.90950601 1.62687398 1.43875694 0.03109259 -2.37265779 -3.45057765
0.19965851 1.72320248 2.2898178 1.34013775 0.40188506 -2.25368524
-4.79999944 -0.29638363 1.49036912 2.26290378 1.88161992 0.70066741
-2.00055299 -4.2727996 -0.37954819 2.06048067 3.136022 2.91270888
0.68341221 -2.54241374 -4.11982953 -0.25818249 2.16914822 3.46019775
2.91502288 1.04311452 -2.62395441 -3.82618485 0.00951645 2.95514821
3.92300846 3.37200094 1.280833 -3.34939407]
332.571174789
我将不胜感激任何解释或猜测默认(Levenberg-Marquardt 方法)可能失败的原因,或者是否有办法对其进行调整以克服这个问题。