scipy 基于其他参数的参数约束非线性回归

nkcskrwz  于 2022-11-10  发布在  其他
关注(0)|答案(1)|浏览(122)

我试图拟合以下形式的非线性回归(两阶段指数衰减):
跨度快速=(Y 0-平台期)* 快速百分比 * .01
慢速跨度=(Y 0-平台期)(100-快速百分比) .01
Y=平台期+跨度快速 * 预期(-KFast * X)+跨度慢速 * 预期(-KSlow * X)
当我在python中使用scipy.optimize.curve_fit进行计算时,我得到的结果与graphpad prism不同,我认为这是因为我没有约束条件:KFast〉KSlow〉0。但是在python中如何做到这一点呢?如何根据另一个参数KFast〉KSlow来约束一个参数。

import scipy.optimize as op
def phaseDecay(x, Y0, Plateau, PercentFast, KFast, KSlow):
    """
    Exponential two phase decay

    Parameters:
    -----------
    """
    SpanFast = (Y0-Plateau)*PercentFast*0.01
    SpanSlow = (Y0-Plateau)*(100-PercentFast)*0.01
    Y = Plateau+SpanFast*np.exp(-KFast*x)+SpanSlow*np.exp(-KSlow*x)
    return Y

popt, pcov = op.curve_fit(
    phaseDecay, 
    df["x"], 
    df["y"],
)

给出:

PARAMETERS:
Y0 100.000000000216
Plateau 69.27241846348228
PercentFast 1.0
KFast 1.0
KSlow 1.0

然而棱镜给出

Y0 100.0
     Plateau    63.58
     PercentFast    72.23
     KFast  0.001626
     KSlow  0.0001125

数据来源:

0       0  100.000000
1    1320   75.323025
2    4500   71.880123
3    7800   70.038842
4   18660   66.408841
5       0  100.000000
6    1500   73.127293
7    4140   68.821849
8    7320   65.775435
9   18540   62.800071
10      0  100.000000
11   1740   75.241496
12   3960   68.779365
13   7440   67.843209
14  18360   65.229471
ttisahbt

ttisahbt1#

如果您将目标函数转换为差的平方和,然后使用支持线性不等式的方法(如SLSQP),则可以使用SciPy来实现这一点。
我是这样做的:

import scipy.optimize as op
import numpy as np

def phaseDecay(x, xdata, ydata):
    """
    Exponential two phase decay

    Parameters:
    -----------
    """

    Y0, Plateau, PercentFast, KFast, KSlow = x

    SpanFast = (Y0-Plateau)*PercentFast*0.01
    SpanSlow = (Y0-Plateau)*(100-PercentFast)*0.01
    Y = Plateau+SpanFast*np.exp(-KFast*xdata)+SpanSlow*np.exp(-KSlow*xdata)
    return np.sum((Y - ydata)**2.0)

def constraints(x):

    return x[-2] - x[-1]

data = np.loadtxt('data_phase_decay.txt')
X0, Y0 = data[:, 1:].T

lb = np.zeros((5, )) + 1e-5
ub = np.array([100.0, 100.0, 100.0, 0.01, 0.01])
x0 = (ub + lb)/2.0
bounds = zip(lb, ub)

# Tell SLSQP that KSlow < KFast

cons = ({'type': 'ineq', 'fun': constraints})

res = op.minimize(phaseDecay, x0, method='SLSQP', constraints=cons,
                  args=(X0, Y0))
print(res)

我得到了以下结果:

fun    : 30.45577816182897
jac    : array([-2.44140625e-04,   1.53160095e-03,  -2.43663788e-04,
                0.00000000e+00,  -2.11278069e+01,   0.00000000e+00])
message: 'Optimization terminated successfully.'
nfev   : 247
nit    : 27
njev   : 27
status : 0
success: True
x      : array([9.99999153e+01,   6.46585138e+01,   6.27658587e+01,
                2.33870123e-02,   2.02987854e-04])

Y0          = 99.999915
Plateau     = 64.658514
PercentFast = 62.765859
KFast       = 0.023387
KSlow       = 0.000203

相关问题