我必须使用 r^2 值将 sigmoid 函数拟合到两种类型的数据集（每种数据集 6 个数据集），然后绘制拟合函数的导数。我收到错误

Question

该代码适用于其中一个数据集，但不适用于另一个数据集。我收到错误（显示在我提供的代码的底部）。

您可以看到我按类型有两个注释掉的数据部分：一个用于工作密度，一个用于应变。 我正在尝试将每种类型中存在的 6 个数据集绘制为针对 x 数据的独立曲线。

我不需要绘制所有 12 个数据 - 我只是根据需要将应变数据交换为工作密度数据。

我无法成功。我确信代码是正确的并且应该可以工作，但是有些东西计算量太大......？

预期结果：6 个数据集的拟合函数（逻辑 sigmoid）图，然后是包含其理论曲线导数的子图。返回 R^2 平方值的值（我不能使用均方误差）。

# This code produces a 99-100% fitness of fit to Normal Fibril strain v. temperature data, or, work      
density v. temperature data, using a technical 
# 4-parameter logistic sigmoid function, but with the centroid parameter pre-defined at 77.5 degrees  
Celsius on the x-axis. This 'b' parameter is 
# left undefined for the derivative function. 

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

# 1. Normal Fibrils user dataset: 
x_data = np.array([25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115,      
120, 125, 130]) 
y_data50 = np.array([0, 0, 0.250155, 1.00062, 1.751085, 2.50155, 3.252015, 4.752945, 6.253875, 
7.00434, 8.50527, 10.756665, 13.00806, 16.260075, 18.51147, 19.261935, 19.51209, 19.762245, 19.762245, 
20.0124, 20.0124, 20.0124]) 
y_data100 = np.array([0, 0, 0, 0, 1.98162, 3.467835, 5.94486, 7.431075, 9.9081, 13.87134, 16.348365,  
18.82539, 22.78863, 28.73349, 34.67835, 40.62321, 42.60483, 43.59564, 44.58645, 45.081855, 45.081855, 
45.081855])
y_data150 = np.array([0, 0, 0, 0, 2.221965, 4.44393, 7.40655, 10.36917, 13.33179, 17.035065, 22.21965,  
25.18227, 31.10751, 37.03275, 44.4393, 51.84585, 62.21502, 68.14026, 71.10288, 74.0655, 75.54681, 
75.54681])
y_data200 = np.array([0, 0, 0, -1.97181, -1.97181, 3.94362, 7.88724, 13.80267, 18.732195, 25.63353,  
31.54896, 39.4362, 45.35163, 51.26706, 62.112015, 73.942875, 86.75964, 96.61869, 106.47774, 108.44955, 
110.42136, 112.39317])
y_data250 = np.array([0, 0, 0, 2.46231, 4.92462, 9.84924, 14.77386, 20.929635, 27.08541, 34.47234, 
41.85927, 54.17082, 64.02006, 73.8693, 86.18085, 103.41702, 118.19088, 130.50243, 137.88936,  
140.35167, 142.81398, 145.27629])
y_data300 = np.array([0, 0, 1.476405, 2.95281, 4.429215, 7.382025, 13.287645, 20.66967, 28.051695,  
38.38653, 51.674175, 67.91463, 82.67868, 100.39554, 121.06521, 137.305665, 150.59331, 157.975335,  
165.35736, 168.31017, 171.26298, 171.26298])

# Normal Fibrils' Work Density Data: 
# m(w) = 50m(f): y_data50 = np.array([0, 0, 0.250155, 1.00062, 1.751085, 2.50155, 3.252015, 4.752945,   
6.253875, 7.00434, 8.50527, 10.756665, 13.00806, 16.260075, 18.51147, 19.261935, 19.51209, 19.762245, 
19.762245, 20.0124, 20.0124, 20.0124]) 
# m(w) = 100m(f): y_data100 = np.array([0, 0, 0, 0, 1.98162, 3.467835, 5.94486, 7.431075, 9.9081, 
13.87134, 16.348365, 18.82539, 22.78863, 28.73349, 34.67835, 40.62321, 42.60483, 43.59564, 44.58645, 
45.081855, 45.081855, 45.081855])
# m(w) = 150m(f): y_data150 = np.array([0, 0, 0, 0, 2.221965, 4.44393, 7.40655, 10.36917, 13.33179,  
17.035065, 22.21965, 25.18227, 31.10751, 37.03275, 44.4393, 51.84585, 62.21502, 68.14026, 71.10288, 
74.0655, 75.54681, 75.54681])
# m(w) = 200m(f): y_data200 = np.array([0, 0, 0, -1.97181, -1.97181, 3.94362, 7.88724, 13.80267, 
18.732195, 25.63353, 31.54896, 39.4362, 45.35163, 51.26706, 62.112015, 73.942875, 86.75964, 96.61869, 
106.47774, 108.44955, 110.42136, 112.39317])
# m(w) = 250m(f): y_data250 = np.array([0, 0, 0, 2.46231, 4.92462, 9.84924, 14.77386, 20.929635, 
27.08541, 34.47234, 41.85927, 54.17082, 64.02006, 73.8693, 86.18085, 103.41702, 118.19088, 130.50243, 
137.88936, 140.35167, 142.81398, 145.27629])
# m(w) = 300m(f): y_data300 = np.array([0, 0, 1.476405, 2.95281, 4.429215, 7.382025, 13.287645,  
20.66967, 28.051695, 38.38653, 51.674175, 67.91463, 82.67868, 100.39554, 121.06521, 137.305665, 
150.59331, 157.975335, 165.35736, 168.31017, 171.26298, 171.26298])

# Normal Fibrils' Strain Data: 
# x_data = np.array([25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115, 
120, 125, 130])
# m(w) = 50m(f): y_data50 = np.array([0, 0, 0.00617284, 0.024691358, 0.043209877, 0.061728395, 
0.080246914, 0.117283951, 0.154320988, 0.172839506, 0.209876543, 0.265432099, 0.320987654, 
0.401234568, 0.456790123, 0.475308642, 0.481481481, 0.487654321, 0.487654321, 0.49382716, 0.49382716, 
0.49382716])
# m(w) = 100m(f): y_data100 = np.array([0, 0, 0, 0, 0, 0, 0.021978022, 0.038461538, 0.065934066, 
0.082417582, 0.10989011, 0.153846154, 0.181318681, 0.208791209, 0.252747253, 0.318681319, 0.384615385, 
0.450549451, 0.472527473, 0.483516484, 0.494505495, 0.5])
# m(w) = 150m(f): y_data150 = np.array([0, 0, 0, 0, 0.015151515, 0.03030303, 0.050505051, 0.070707071, 
0.090909091, 0.116161616, 0.151515152, 0.171717172, 0.212121212, 0.252525253, 0.303030303, 
0.353535354, 0.424242424, 0.464646465, 0.484848485, 0.505050505, 0.515151515, 0.515151515])
# m(w) = 200m(f): y_data200 = np.array([0, 0, 0, -0.009345794, -0.009345794, 0.018691589, 0.037383178,     
0.065420561, 0.088785047, 0.121495327, 0.14953271, 0.186915888, 0.214953271, 0.242990654, 0.294392523,  
0.35046729, 0.411214953, 0.457943925, 0.504672897, 0.514018692, 0.523364486, 0.53271028])
# m(w) = 250m(f): y_data250 = np.array([0, 0, 0, 0.009090909, 0.018181818, 0.036363636, 0.054545455,  
0.077272727, 0.1, 0.127272727, 0.154545455, 0.2, 0.236363636, 0.272727273, 0.318181818, 0.381818182, 
0.436363636, 0.481818182, 0.509090909, 0.518181818, 0.527272727, 0.536363636])
# m(w) = 300m(f): y_data300 = np.array([0, 0, 0.004504505, 0.009009009, 0.013513514, 0.022522523, 
0.040540541, 0.063063063, 0.085585586, 0.117117117, 0.157657658, 0.207207207, 0.252252252, 
0.306306306, 0.369369369, 0.418918919, 0.459459459, 0.481981982, 0.504504505, 0.513513514, 
0.522522523, 0.522522523])

# 2. Four-parameter sigmoid fitting function: 
def sigmoid(x, a, b, c, U):
    return (U / (1 + np.exp(-c * (x - b))))**a

# 3. Sigmoid derivative function: 
def sigmoid_derivative(x, a, b, c, U):
    return (a * U * np.exp(-c * (x - b)) * (1 + np.exp(-c * (x - b)))**(-a - 1) * c) / (1 + np.exp(-c  
* (x - b)))**(2*a)

# 4. Initial guess for the fitting parameters: 
initial_guess = [1, 77.5, 0.1, 0.5]

# 5. Fit the sigmoid function to each dataset and calculate R-squared value: 
popt50, _ = curve_fit(sigmoid, x_data, y_data50, p0=initial_guess)
y_pred50 = sigmoid(x_data, *popt50)
r_squared_50 = r2_score(y_data50, y_pred50)

popt100, _ = curve_fit(sigmoid, x_data, y_data100, p0=initial_guess)
y_pred100 = sigmoid(x_data, *popt100)
r_squared_100 = r2_score(y_data100, y_pred100)

popt150, _ = curve_fit(sigmoid, x_data, y_data150, p0=initial_guess)
y_pred150 = sigmoid(x_data, *popt150)
r_squared_150 = r2_score(y_data150, y_pred150)

popt200, _ = curve_fit(sigmoid, x_data, y_data200, p0=initial_guess)
y_pred200 = sigmoid(x_data, *popt200)
r_squared_200 = r2_score(y_data200, y_pred200)

popt250, _ = curve_fit(sigmoid, x_data, y_data250, p0=initial_guess)
y_pred250 = sigmoid(x_data, *popt250)
r_squared_250 = r2_score(y_data250, y_pred250)

popt300, _ = curve_fit(sigmoid, x_data, y_data300, p0=initial_guess)
y_pred300 = sigmoid(x_data, *popt300)
r_squared_300 = r2_score(y_data300, y_pred300)

# 6. Generate x values for the fitted curve: 
x_fit = np.linspace(min(x_data), max(x_data), 1000)

# 7. Calculate the y values using the fitted parameters for each dataset: 
y_fit50 = sigmoid(x_fit, *popt50)
y_fit100 = sigmoid(x_fit, *popt100)
y_fit150 = sigmoid(x_fit, *popt150)
y_fit200 = sigmoid(x_fit, *popt200)
y_fit250 = sigmoid(x_fit, *popt250)
y_fit300 = sigmoid(x_fit, *popt300)

# 8. Calculate the derivatives of each fitted curve: 
dy_dx_fit50 = sigmoid_derivative(x_fit, *popt50)
dy_dx_fit100 = sigmoid_derivative(x_fit, *popt100)
dy_dx_fit150 = sigmoid_derivative(x_fit, *popt150)
dy_dx_fit200 = sigmoid_derivative(x_fit, *popt200)
dy_dx_fit250 = sigmoid_derivative(x_fit, *popt250)
dy_dx_fit300 = sigmoid_derivative(x_fit, *popt300)

# 9. Plotting: 
plt.figure(figsize=(15, 5))

# 10. Plot sigmoid fits: 
plt.subplot(1, 2, 1) 
plt.plot(x_fit, y_fit300, '-', color = 'black', label = f'm(w) = 300m(f)|(R^2={r_squared_300:.2f})') 
plt.plot(x_fit, y_fit250, '-', color = 'purple', label = f'm(w) = 250m(f)|(R^2={r_squared_250:.2f})')
plt.plot(x_fit, y_fit200, '-', color = 'red', label = f'm(w) = 200m(f)|(R^2={r_squared_200:.2f})')
plt.plot(x_fit, y_fit150, '-', color = 'orange', label = f'm(w) = 150m(f)|(R^2={r_squared_150:.2f})')
plt.plot(x_fit, y_fit100, '-', color = 'gold', label = f'm(w) = 100m(f)|(R^2={r_squared_100:.2f})')
plt.plot(x_fit, y_fit50, '-', color = 'green', label = f'm(w) = 50m(f)|(R^2={r_squared_50:.2f})')

plt.xlabel('Temperature (°C)') 
plt.ylabel('Work Density WD') 
plt.title('Sigmoidal Fit to Normal Fibril Data [WD v. T]') 
#plt.ylabel('Strain ε, (Δl/$\mathregular{l_{0}}$)') 
#plt.title('Sigmoidal Fit to Normal Fibril Data [ε v. T]') 
plt.legend() 
plt.grid(True)

# 11. Plot derivatives: 
plt.subplot(1, 2, 2)
plt.plot(x_fit, dy_dx_fit300, '--', color = 'black', label = 'Rate for m(w) = 300m(f)')
plt.plot(x_fit, dy_dx_fit250, '--', color = 'purple', label = 'Rate for m(w) = 250m(f)')
plt.plot(x_fit, dy_dx_fit200, '--', color = 'red', label = 'Rate for m(w) = 200m(f)')
plt.plot(x_fit, dy_dx_fit150, '--', color = 'orange', label = 'Rate for m(w) = 150m(f)')
plt.plot(x_fit, dy_dx_fit100, '--', color = 'gold', label = 'Rate for m(w) = 100m(f)')
plt.plot(x_fit, dy_dx_fit50, '--', color = 'green', label='Rate for m(w) = 50m(f)')

plt.xlabel('Temperature (°C)') 
plt.ylabel('ΔWD/ΔT, (Nm/kg°C)')
#plt.ylabel('Δε/ΔT, ((Δl/$\mathregular{l_{0}}$)/°C)')
plt.title('Contraction Rate')
plt.legend()
plt.grid(True)
plt.show() 

错误：

“运行时警告：幂中遇到除零”
“运行时警告：电源中遇到无效值”
“优化警告：无法估计参数的协方差。”
'第 52 行，在 popt100，_ ..... 第 982 行，在 curve_fit 中引发运行时错误（“未找到最佳参数:)”
“未找到最佳参数：函数调用次数已达到 maxfev = 1000。”

Answer 1

您面临的问题来自于 sigmoid 函数试图将负数乘以分数次方，这会导致 Python 中的复数。我建议调整参数的界限，以确保输入幂函数的所有值都是正数，修改参数的初始猜测，以确保它们更有可能产生有效的解决方案，并添加界限来拟合曲线。

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


x_data = np.array([25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90, 95, 100, 105, 110, 115, 120, 125, 130])

y_datasets = {
    '50': np.array([0, 0, 0.250155, 1.00062, 1.751085, 2.50155, 3.252015, 4.752945, 6.253875, 
                    7.00434, 8.50527, 10.756665, 13.00806, 16.260075, 18.51147, 19.261935, 19.51209, 19.762245, 19.762245, 
                    20.0124, 20.0124, 20.0124]), 
    '100': np.array([0, 0, 0, 0, 1.98162, 3.467835, 5.94486, 7.431075, 9.9081, 13.87134, 16.348365,  
                     18.82539, 22.78863, 28.73349, 34.67835, 40.62321, 42.60483, 43.59564, 44.58645, 45.081855, 45.081855, 
                     45.081855]),
    '150': np.array([0, 0, 0, 0, 2.221965, 4.44393, 7.40655, 10.36917, 13.33179, 17.035065, 22.21965,  
                     25.18227, 31.10751, 37.03275, 44.4393, 51.84585, 62.21502, 68.14026, 71.10288, 74.0655, 75.54681, 
                     75.54681]),
    '200': np.array([0, 0, 0, -1.97181, -1.97181, 3.94362, 7.88724, 13.80267, 18.732195, 25.63353,  
                     31.54896, 39.4362, 45.35163, 51.26706, 62.112015, 73.942875, 86.75964, 96.61869, 106.47774, 108.44955, 
                     110.42136, 112.39317]),
    '250': np.array([0, 0, 0, 2.46231, 4.92462, 9.84924, 14.77386, 20.929635, 27.08541, 34.47234, 
                     41.85927, 54.17082, 64.02006, 73.8693, 86.18085, 103.41702, 118.19088, 130.50243, 137.88936,  
                     140.35167, 142.81398, 145.27629]),
    '300': np.array([0, 0, 1.476405, 2.95281, 4.429215, 7.382025, 13.287645, 20.66967, 28.051695,  
                     38.38653, 51.674175, 67.91463, 82.67868, 100.39554, 121.06521, 137.305665, 150.59331, 157.975335,  
                     165.35736, 168.31017, 171.26298, 171.26298])
}

def sigmoid(x, a, c, U):
    b = 77.5  
    base = U / (1 + np.exp(-c * (x - b)))
    return np.maximum(base, 0) ** a  

def sigmoid_derivative(x, a, c, U):
    b = 77.5
    e = np.exp(-c * (x - b))
    base = U / (1 + e)
    return (a * c * e * base ** (a - 1)) / (1 + e) ** a

bounds = ([0, 0, 0], [np.inf, 1, np.inf])  
fit_results = {}
for key, y_data in y_datasets.items():
    try:
        popt, _ = curve_fit(sigmoid, x_data, y_data, p0=[1, 0.1, max(y_data)], bounds=bounds, maxfev=5000)
        y_pred = sigmoid(x_data, *popt)
        r_squared = r2_score(y_data, y_pred)
        fit_results[key] = (popt, r_squared)
    except Exception as e:
        print(f"Error fitting dataset {key}: {e}")

x_fit = np.linspace(min(x_data), max(x_data), 1000)

plt.figure(figsize=(15, 10))

plt.subplot(2, 1, 1)
colors = ['green', 'gold', 'orange', 'red', 'purple', 'black']
for i, (key, (popt, r_squared)) in enumerate(fit_results.items()):
    y_fit = sigmoid(x_fit, *popt)
    plt.plot(x_fit, y_fit, '-', color=colors[i], label=f'm(w) = {key}m(f) | R^2={r_squared:.2f}')

plt.xlabel('Temperature (°C)')
plt.ylabel('Work Density or Strain')
plt.title('Sigmoidal Fit to Data')
plt.legend()
plt.grid(True)

plt.subplot(2, 1, 2)
for i, (key, (popt, _)) in enumerate(fit_results.items()):
    dy_dx_fit = sigmoid_derivative(x_fit, *popt)
    plt.plot(x_fit, dy_dx_fit, '--', color=colors[i], label=f'Rate for m(w) = {key}m(f)')

plt.xlabel('Temperature (°C)')
plt.ylabel('Rate of Change')
plt.title('Derivative of Sigmoidal Fit')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

这解决了问题并给出了