如何在MATLAB中线性化数据并寻找最优b使log(y)与log((x-b)/b)呈线性

Got it, let's tackle this problem step by step. You want to find the optimal parameter b such that log(y) has a linear relationship with log((x - b)/b). Here's how you can do this using Python with numerical optimization and linear regression:

1. 问题数学转化 #

首先，我们的目标是让：
$$\log(y) = k \cdot \log\left(\frac{x - b}{b}\right) + c$$
其中$k$是斜率，$c$是截距。等价于非线性模型：
$$y = A \cdot \left(\frac{x - b}{b}\right)^k$$
（其中$A = e^c$）

我们需要找到最小化$\log(y)$与拟合值之间残差平方和的$b$。

2. 代码实现 #

导入必要库 #

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

定义原始数据 #

x = np.array([45.5976, 45.6311, 45.6599, 45.679, 45.703, 45.7461, 45.7749])
y = np.array([0.17, 1.7, 5.1, 17, 51, 170, 510])

定义目标优化函数 #

这个函数会计算给定b时，$\log(y)$与$\log((x-b)/b)$线性拟合的残差平方和，我们要最小化这个值：

def objective(b):
    # 确保x > b，避免对数运算出现无效值（负数/零）
    if np.any(x <= b):
        return np.inf  # 返回无穷大，让优化器自动避开这种情况
    
    # 计算转换后的自变量和因变量
    log_transformed_x = np.log((x - b) / b)
    log_y = np.log(y)
    
    # 拟合线性模型，得到斜率和截距
    slope, intercept = np.polyfit(log_transformed_x, log_y, 1)
    
    # 计算预测值与真实值的残差平方和
    predicted_log_y = slope * log_transformed_x + intercept
    return np.sum((log_y - predicted_log_y) ** 2)

寻找最优b #

我们需要给b一个初始猜测值（必须小于最小的x值，这里最小x是45.5976，所以选45.5作为初始值）：

# 初始猜测值
initial_b = 45.5

# 使用Nelder-Mead优化算法（无需导数，适合这类问题）
optimization_result = minimize(objective, initial_b, method='Nelder-Mead')

# 提取最优b
optimal_b = optimization_result.x[0]
print(f"找到的最优参数b: {optimal_b:.6f}")
# 可以检查优化是否成功
print(f"优化是否收敛: {optimization_result.success}")

验证拟合效果 #

# 用最优b计算转换后的变量
log_transformed_x_opt = np.log((x - optimal_b) / optimal_b)
log_y = np.log(y)

# 重新拟合线性模型
slope_opt, intercept_opt = np.polyfit(log_transformed_x_opt, log_y, 1)
predicted_log_y_opt = slope_opt * log_transformed_x_opt + intercept_opt

# 绘制对比图
plt.figure(figsize=(12, 6))

# 原始数据非线性曲线
plt.subplot(1, 2, 1)
plt.plot(x, y, 'o', color='tab:blue', label='原始数据点')
plt.xlabel('x')
plt.ylabel('y')
plt.title('原始数据的非线性关系')
plt.legend()

# 转换后的线性拟合图
plt.subplot(1, 2, 2)
plt.plot(log_transformed_x_opt, log_y, 'o', color='tab:orange', label='转换后数据点')
plt.plot(log_transformed_x_opt, predicted_log_y_opt, '-', color='tab:green', 
         label=f'拟合直线: log(y) = {slope_opt:.2f}*log((x-b)/b) + {intercept_opt:.2f}')
plt.xlabel('log((x - b)/b)')
plt.ylabel('log(y)')
plt.title(f'最优b={optimal_b:.6f}时的线性拟合')
plt.legend()

plt.tight_layout()
plt.show()

3. 关键注意事项 #

• 初始值选择：b必须严格小于所有x值，否则$\log((x-b)/b)$会出现无意义的负数/零，所以初始猜测要小于最小的x。
• 优化算法：这里用Nelder-Mead是因为它不需要计算目标函数的导数，适合这种无法直接求导的复杂目标。
• 结果验证：一定要检查optimization_result.success是否为True，确保优化收敛。

内容的提问来源于stack exchange，提问作者Nini