SciPy.optimize.least_squares() 5PL曲线优化问题
SciPy.optimize.least_squares() 5PL Curve Optimization problems
我正在尝试编写一个脚本,它将采用 x 和 y 值的输入数组并将它们拟合到 5-PL 曲线(由方程 F(x) = D+(A-D)/((1 +(x/C)^B)^E))。然后我希望能够使用预测曲线来获取给定的 y 值并从曲线中推断出 x 值,由方程 F(y) = C(((A-D)/(-D+y))^ (1/E)-1)^(1/B).
下面的回答修复了之前的错误,但是还是不太合适。我已经引入了一个打印函数,其中包含输入 curve_fit 范围内的少量 y 值,它在整个范围内产生几乎完全相同的 x 值。有什么想法吗?
编辑:对于现在正在看的任何人来说,问题似乎出在我对 B 的估计上。在大多数情况下,山坡应该在 -1 和 1 之间,而不是数千。这使得估计太远了。
import numpy as np
import scipy.optimize as sp
def logistic5(x, A, B, C, D, E):
'''5PL logistic equation'''
log = D + (A-D)/(np.power((1 + np.power((x/C), B)), E))
return log
def residuals(p, y, x):
'''Deviations of data from fitted 5PL curve'''
A, B, C, D, E = p
err = y - logistic5(x, A, B, C, D, E)
print(err)
return err
def log_solve_for_x(curve, y):
'''Returns the estimated x value for the provided y value'''
A, B, C, D, E = curve
return C*(np.power((np.power(((A-D)/(-D+y)), (1/E))-1), (1/B)))
# Toy data set
x = np.array([130, 38, 15, 4.63, 1.41])
y = np.array([9121, 1987, 1017, 343, 117])
# Set initial guess for parameters
A = np.amin(y) # Min asymptote
D = np.amax(y) # Max asymptote
B = (D-A)/(np.amax(x)-np.amin(x)) # Steepness
C = (np.amax(x)-np.amin(x))/2 # inflection point
E = 1 # Asymmetry factor
# Optimize curve for initial parameters
p0 = [A, B, C, D, E]
# set bounds for each parameter
pu = []
pl = []
for p in p0:
pu.append(p*1.5)
pl.append(p*0.5)
print(pu)
print(pl)
print("Initial guess of parameters is: ", p0)
curve = sp.least_squares(fun=residuals, x0=p0, args=(y, x), bounds=(pl, pu))
curve = curve.x.tolist()
print("Optimized curve parameters are: ", curve)
# Predict x values based on given y
y = [1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000]
for sample in y:
solve = log_solve_for_x(curve, sample)
print("Predicted X value for y =", sample, " is: ", solve)
您的曲线未针对任何参数值定义。但是您没有为 least_squares
提供该信息。在某些时候,求解器进入一个不允许的区域并卡在那里从 residuals
获取 nans,你会收到有关无效功率的消息。您拥有微不足道的权力,可能只是设置 E>=0, B>=0
。但是你的基础是不平凡的。您需要切换到支持通用约束(例如 scipy.optimize.minimize
)的求解器并添加 base >=0
的约束,或者以其他方式将搜索限制在可接受的域中,例如:
pu = []
pl = []
for p in p0:
pu.append(p*1.5)
pl.append(p*.5)
curve = sp.least_squares(fun=residuals, x0=p0, args=(y, x), bounds=(pl, pu))
您也可以尝试修复残差,使其适用于任何参数,例如用到初始猜测的距离替换 nan。但它可能工作效率低下。
要改善拟合结果,您可以尝试更好的初始点或多起点或两者。
A = np.amin(y) # Min asymptote
D = np.amax(y) # Max asymptote
B = (D-A)/np.amax(x)*10 # Steepness
C = np.amax(x)/10 # inflection point
E = 0.001 # Asymmetry factor
p0 = [A, B, C, D, E]
print("Initial guess of parameters is: ", p0)
pu = []
pl = []
for p in p0:
pu.append(p*1.5)
pl.append(p*.5)
best_cost = np.inf
for i in range(100):
for i in range(5):
p0[i] = np.random.uniform(pl[i], pu[i])
curve = sp.least_squares(fun=residuals, x0=p0, args=(y, x), bounds=(pl, pu))
print(p0, curve.cost)
if best_cost > curve.cost:
best_cost = curve.cost
curve_out = curve.x.tolist()
print("Optimized curve parameters are: ", curve_out)
plt.plot(x, y, '.')
xx = np.linspace(0, 150, 100)
yy = []
for x in xx:
yy.append(logistic5(x, *curve_out))
plt.plot(xx, yy)
plt.show()
我正在尝试编写一个脚本,它将采用 x 和 y 值的输入数组并将它们拟合到 5-PL 曲线(由方程 F(x) = D+(A-D)/((1 +(x/C)^B)^E))。然后我希望能够使用预测曲线来获取给定的 y 值并从曲线中推断出 x 值,由方程 F(y) = C(((A-D)/(-D+y))^ (1/E)-1)^(1/B).
下面的回答修复了之前的错误,但是还是不太合适。我已经引入了一个打印函数,其中包含输入 curve_fit 范围内的少量 y 值,它在整个范围内产生几乎完全相同的 x 值。有什么想法吗?
编辑:对于现在正在看的任何人来说,问题似乎出在我对 B 的估计上。在大多数情况下,山坡应该在 -1 和 1 之间,而不是数千。这使得估计太远了。
import numpy as np
import scipy.optimize as sp
def logistic5(x, A, B, C, D, E):
'''5PL logistic equation'''
log = D + (A-D)/(np.power((1 + np.power((x/C), B)), E))
return log
def residuals(p, y, x):
'''Deviations of data from fitted 5PL curve'''
A, B, C, D, E = p
err = y - logistic5(x, A, B, C, D, E)
print(err)
return err
def log_solve_for_x(curve, y):
'''Returns the estimated x value for the provided y value'''
A, B, C, D, E = curve
return C*(np.power((np.power(((A-D)/(-D+y)), (1/E))-1), (1/B)))
# Toy data set
x = np.array([130, 38, 15, 4.63, 1.41])
y = np.array([9121, 1987, 1017, 343, 117])
# Set initial guess for parameters
A = np.amin(y) # Min asymptote
D = np.amax(y) # Max asymptote
B = (D-A)/(np.amax(x)-np.amin(x)) # Steepness
C = (np.amax(x)-np.amin(x))/2 # inflection point
E = 1 # Asymmetry factor
# Optimize curve for initial parameters
p0 = [A, B, C, D, E]
# set bounds for each parameter
pu = []
pl = []
for p in p0:
pu.append(p*1.5)
pl.append(p*0.5)
print(pu)
print(pl)
print("Initial guess of parameters is: ", p0)
curve = sp.least_squares(fun=residuals, x0=p0, args=(y, x), bounds=(pl, pu))
curve = curve.x.tolist()
print("Optimized curve parameters are: ", curve)
# Predict x values based on given y
y = [1000, 2000, 3000, 4000, 5000, 6000, 7000, 8000, 9000]
for sample in y:
solve = log_solve_for_x(curve, sample)
print("Predicted X value for y =", sample, " is: ", solve)
您的曲线未针对任何参数值定义。但是您没有为 least_squares
提供该信息。在某些时候,求解器进入一个不允许的区域并卡在那里从 residuals
获取 nans,你会收到有关无效功率的消息。您拥有微不足道的权力,可能只是设置 E>=0, B>=0
。但是你的基础是不平凡的。您需要切换到支持通用约束(例如 scipy.optimize.minimize
)的求解器并添加 base >=0
的约束,或者以其他方式将搜索限制在可接受的域中,例如:
pu = []
pl = []
for p in p0:
pu.append(p*1.5)
pl.append(p*.5)
curve = sp.least_squares(fun=residuals, x0=p0, args=(y, x), bounds=(pl, pu))
您也可以尝试修复残差,使其适用于任何参数,例如用到初始猜测的距离替换 nan。但它可能工作效率低下。
要改善拟合结果,您可以尝试更好的初始点或多起点或两者。
A = np.amin(y) # Min asymptote
D = np.amax(y) # Max asymptote
B = (D-A)/np.amax(x)*10 # Steepness
C = np.amax(x)/10 # inflection point
E = 0.001 # Asymmetry factor
p0 = [A, B, C, D, E]
print("Initial guess of parameters is: ", p0)
pu = []
pl = []
for p in p0:
pu.append(p*1.5)
pl.append(p*.5)
best_cost = np.inf
for i in range(100):
for i in range(5):
p0[i] = np.random.uniform(pl[i], pu[i])
curve = sp.least_squares(fun=residuals, x0=p0, args=(y, x), bounds=(pl, pu))
print(p0, curve.cost)
if best_cost > curve.cost:
best_cost = curve.cost
curve_out = curve.x.tolist()
print("Optimized curve parameters are: ", curve_out)
plt.plot(x, y, '.')
xx = np.linspace(0, 150, 100)
yy = []
for x in xx:
yy.append(logistic5(x, *curve_out))
plt.plot(xx, yy)
plt.show()