2013-06-24 34 views
3
numseq = ['0012000', '0112000', '0212000', '0312000', '1012000', '1112000',                     '1212000', '1312000', '2012000', '2112000', '2212000', '2312000', '3012000', '3112000',   '3212000', '3312000', '0002000', '0022000', '0032000', '1002000', '1022000', '1032000',  '2002000', '2022000', '2032000', '3002000', '3022000', '3032000', '0010000', '0011000', '0013000', '1010000', '1011000', '1013000', '2010000', '2011000', '2013000', '3010000', '3011000', '3013000', '0012100', '0012200', '0', '1012100', '1012200', '1', '2012100', '2012200', '2', '3012100'] 
prob = [-0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.78361598908750163, -0.66474525640568083, -0.66474525640568083, -0.66474525640568083, -0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.66474525640568083, -0.66474525640568083, -0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.66474525640568083, -0.66474525640568083, -0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.66474525640568083, -0.66474525640568083, -0.66474525640568083, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212, -0.49518440694747212] 

numseqprob是每個長度爲50的列表。他們是收集的實驗數據。 numseq對應於X軸值,並且prob對應於Y軸值。scipy optimize fmin語法

,我希望儘量減少的功能是:

def residue(allparams, xdata, ydata): 
    chi2 = 0.0 
    for i in range(0,len(xdata)): 
     x = xdata[i] 
     y = 0 
     for j in range(len(x)): 
      y = y-allparams[int(x[j])][j] 
      chi2 = chi2 + (ydata[i]-y)**2 
return chi2 

所以:

  • allparams是一個4×7的矩陣,其中包含的所有參數進行優化。
  • xdata是X軸值,即numseq
  • ydata僅僅是一個數字的列表即prob

chi2是實驗和模型值之間的平方差。這是必須最小化的。

x0 = [[-0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6], [-0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6], [-0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6], [-0.6, -0.6, -0.6, -0.6, -0.6, -0.6, -0.6]] 

現在我該怎樣調用這個函數fmin

初始猜測的參數由下式給出?我試圖

fmin(residue, x0, args=(numseq, prob)) 

,但我不斷收到一個錯誤:

Traceback (most recent call last): 
    File "<pyshell#362>", line 1, in <module> 
    fmin(residue, x0, args=(numseq, prob)) 
    File "C:\Python31\lib\site-packages\scipy\optimize\optimize.py", line 258, in fmin 
    fsim[0] = func(x0) 
    File "C:\Python31\lib\site-packages\scipy\optimize\optimize.py", line 177, in function_wrapper 
    return function(x, *args) 
    File "<pyshell#361>", line 7, in residue 
    y = y-allparams[int(x[j])][j] 
IndexError: invalid index to scalar variable. 

爲什麼會這樣?是否因爲fmin不能接受二維數組作爲初始猜測?那麼我是否必須將我的整個代碼更改爲一維數組參數?

即使您無法解釋此問題,您是否至少可以告訴我fmin模塊的工作原理?即如何實現fmin用於優化N維數組的語法?你能解釋一下args()是什麼嗎?我是新來的優化,我不知道如何實現它:(

+0

首先發布的'residue'功能似乎不完整。 – Simon

回答

2

「fmin」例程可以接受2d數組作爲初始猜測。但它做的第一件事是平坦這個數組[[4,所以會發生的是,你的餘數函數需要一個(4,7)數組作爲輸入,而「fmin」例程給它一個長度爲28的扁平「x0」。這就是爲什麼你看到錯誤:。
y = y-allparams[int(x[j])][j]
IndexError: invalid index to scalar variable.

See the source code here.

這樣看來,你將不得不不過更改殘留函數接受一個向量而不是數組這似乎並不太壞我嘗試以下哪似乎工作(注意:請仔細檢查!)

def residue_alternative(allparams, inshape, xdata, ydata): 
    m, n = inshape 
    chi2 = 0.0 
    for i in range(0,len(xdata)): 
     x = xdata[i] 
     y = 0 
     for j in range(len(x)): 
      idx = int(x[j]) * n + j #Double check this to 
      y = y-allparams[idx]  #make sure it does what you want 
      chi2 = chi2 + (ydata[i]-y)**2 
    return chi2 

我把它叫做使用:

x0 = -0.6 * np.ones((4,7), dtype=np.double) 
[xopt, fopt, iter, funcalls, warnflag] = \ 
    fmin(residue_alternative, x0, args=(x0.shape, numseq, prob), 
     maxiter = 100000, 
     maxfun = 100000, 
     full_output=True, 
     disp=True) 

並取得了以下成果:

Optimization terminated successfully. 
     Current function value: 7.750523 
     Iterations: 21570 
     Function evaluations: 26076 

>>>xopt 
array([ 0.57669042, -0.21965861, 0.2635061 , -0.08284016, -0.0779489 , 
    -0.10358114, 0.14041582, 0.72469391, -0.43190214, 0.31269757, 
    -0.0338726 , -0.14919739, -2.58314651, 2.74251214, 0.57695759, 
    -0.49574628, 0.1490926 , 0.04912353, 0.02420988, 1.17924051, 
    -7.2147027 , 0.57860843, -0.28386938, 0.2431877 , -0.22674694, 
    -0.58308225, -6.05706775, -2.06350063])  

,你可以重塑一個4X7陣列。試試這個,讓我知道它是否有效/有幫助。