0

我正在使用 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 方法)可能失败的原因,或者是否有办法对其进行调整以克服这个问题。

4

0 回答 0