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', '0012300', '1012100', '1012200', '1012300', '2012100', '2012200', '2012300', '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]
numseq
和 prob
都是长度为 50 的列表。它们是收集的实验数据。 numseq
对应X轴值,prob
对应Y轴值。
我要最小化的函数是:
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()
是什么吗?我是优化的新手,我不知道如何实现它:(
最佳答案
“fmin”例程可以接受二维数组作为初始猜测。但它做的第一件事是展平这个数组 [(4,7) --> (28)]。所以发生的是你的残差函数将 (4,7) 数组作为输入,“fmin”例程给它一个长度为 28 的扁平化“x0”。这就是你看到错误的原因: y = y-allparams[int(x[j])][j]
IndexError:标量变量的索引无效。
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])
您可以将其 reshape 为 4x7 阵列。试一试,让我知道它是否有效/有帮助。
关于python - scipy 优化 fmin 语法,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/17269149/