拟合指数衰减

Question

我正在尝试求解以下线性化方程：

ln⁡{1−y/y}=ln⁡(c)−b(x)

使用 python scipy curvefit 或其他类似方法，请告诉我该怎么做？

样本数据：

x = [15, 16, 17, 18, 19, 20]

y = [0.78, 0.67, 0.56, 0.41, 0.31, 0.20]

到目前为止我尝试过的代码：

import numpy as np
import scipy as sp
import pylab
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import math
import warnings   

def sigmoid(x,c,b):
    y = np.log(c)-b*x
    return y

def sigmoid_solve(y, c, b):
    x = (np.log(c)+np.log((1-y)/y))/b
    return x

y_new = []

x_data = [15, 16, 17, 18, 19, 20]
y_data = [0.78, 0.67, 0.56, 0.41, 0.31, 0.20]


for data in y_data:
    y_new.append(np.log((1-data)/data))


popt, pcov = curve_fit(sigmoid, x_data, y_new)
ce50 = sigmoid_solve(0.5,popt[0],popt[1])

x = np.linspace(10,40,10)
y = 1/(1+np.exp(sigmoid(x, *popt)))
plt.plot(x, y, 'r',label='logistic exp curve fit')
plt.plot(x_data, y_data,'o',label='data plot')
plt.ylim(0, 1)
plt.xlim(10, 50)
plt.legend(loc='best')
plt.savefig('test.png')
plt.close("all")

有没有更好的方法来解决这个问题？

Answer 1

我不会很好地理解你的问题。

在这一刻，我会思考很多方面。 （当然我不知道他们是不是最好的）

1.- 我会接受你的功能。

ln⁡{1−y/y}=ln⁡(c)−b(x)

Anb 如果这是您要求解的方程（找到曲线的交点），我会找到 ln⁡{1−y/y}-(ln⁡(c)−b(x)) 的解=0。

2.-如果你想知道拟合数据的最佳曲线，你会使用贝塞尔曲线的拟合，matplotlib已经有一些功能，看这里：

Bézier example

3.- 但是如果你想解决这个系统（我认为它是一个 EDO）你会在这里检查，python 有很多这样的：

Scipy Integration and ODEs

还有这个：

Solve Differential Equations in Python

干杯。

安东尼奥。

Answer 2

我怀疑您代码中的问题是 c 参数，因为 curve_fit 尝试使用负值导致错误。不幸的是，我在文档中没有找到任何关于如何“绑定”输出参数的参考（有人可以通过指出如何做到这一点来改进这个答案）。尝试用 log_c 替换 c，它可以有任何值：

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

def sigmoid(x, log_c, b):
    return log_c - b * x

def sigmoid_solve(y, log_c, b):
    x = (log_c + np.log((1-y)/y))/b
    return x

x_data = np.array([15, 16, 17, 18, 19, 20])
y_data = np.array([0.78, 0.67, 0.56, 0.41, 0.31, 0.20])
y_new = (np.log((1.0-y_data)/y_data))

popt, pcov = curve_fit(sigmoid, x_data, y_new)
ce50 = sigmoid_solve(0.5,popt[0],popt[1])

x = np.linspace(10,40,10)
y = 1/(1+np.exp(sigmoid(x, *popt)))

plt.plot(x, y, 'r',label='logistic exp curve fit')
plt.plot(x_data, y_data,'o',label='data plot')
plt.ylim(0, 1)
plt.xlim(10, 50)
plt.legend(loc='best')
plt.savefig('test.png')
plt.close("all")

这段代码对我来说很好用。另请注意，如果您使用 np.array 而不是列表，则可以应用不带循环的数学函数。