利用最小二乘法进行数据的精确拟合，如何调节初值使得拟合效果更好

Question

clc;
clear all;
close all;
t = [0 1.5 4.5 10.5 22.5 45 55.5 75 102 135 162 202.5 238.5 277.5 304.5 339 376.5 402 442.5 504 546 600];
y = [0.081571654 0.094432787 0.109133055 0.123423916 0.134230049 0.141051663 0.142588155 0.144447904 0.14618251 0.147662514 0.148371117 0.149273703 0.150370796 0.151081321 0.151469093 0.151986357 0.152568599 0.152827482 0.153345449 0.153993285 0.154252537 0.154900961];
k=0.44;
fun=@(x,t)k/x(1)+k/x(2)*t.^(x (3))./gamma (1+x (3));
%x0 = [100,10,1];
x0 = [80,10,0.1];
x = lsqcurvefit(fun,x0,t,y)
times = linspace(t(1),t(end));
%tempdata=(y-fun(x,times)).^2;
%tempdata2=(y-mean(fun(x,times))).^2;
%r2=1 - ( sum(tempdata)/sum(tempdata2) );
%sqrt(mean((y-fun(x,times)).^2));
%y1为预测值 y为实际值
plot(t,y,'ko',times,fun(x,times),'b-')
legend('Data','Fitted exponential')
xlabel('x')
ylabel('y')
title('Data and Fitted Curve')
R2=1 - (sum((fun(x,t)- y).^2)/ sum((y - mean(y)).^2))
RMSE = sqrt(mean((fun(x,t)-y).^2))

Answer 1

引自免费微信小程序：皆我百晓生

在使用最小二乘法进行数据拟合时，选择合适的初值对于获得更好的拟合效果至关重要。从您提供的代码来看，您正在尝试对数据进行双曲正弦函数（或者是修正的指数函数）类型的拟合，函数形式为 k/x1 + k/x2*t^(x3)/gamma(1+x3)。

针对如何调节初值以优化拟合结果，以下是一些策略：

1. 基于领域知识：如果你对所研究的现象有深入理解，可以根据实际背景设置合理的初值范围。例如，在您的问题中，如果知道衰减率、时间常数或指数项的大致范围，可以据此设定初始猜测值。

2. 观察数据趋势：通过可视化数据并观察其大致形状和变化趋势，可以帮助你估计参数的大致大小。例如，你可以粗略估计一下 t=0 时的截距、斜率以及曲线陡峭程度对应的指数项。

3. 分段拟合：如果数据的变化趋势复杂，可以考虑先对数据的不同部分进行局部拟合，从而获取每个部分参数的合理区间，再综合全局优化。

4. 网格搜索：可以在一个合理的参数空间范围内，采用网格搜索方法遍历不同的初值组合，并记录下每次拟合得到的目标函数（残差平方和）值，选取使目标函数值最小的初值组合。

5. 随机搜索：与网格搜索类似，但不是均匀地在整个搜索空间取点，而是随机生成一组初值，多次迭代后找到最优解。

对于上述代码中的具体例子，您已经提供了 [80,10,0.1] 这样的初值。您可以尝试调整这些值，或者采取上述方法来确定更优的初值。例如，可以通过多次运行程序，每次都更改 x0 的值，然后比较不同初值下的拟合效果（R2分数或RMSE等指标）。

此外，请注意，lsqcurvefit 函数本身具备一定的全局搜索能力，但是它是一个局部优化算法，所以合适的初值会提高找到全局最优解的概率。若数据集较复杂且模型可能存在多个局部极小值，可能需要结合多起始点或多轮优化来确保结果的可靠性。

下面是一个示例，展示如何修改代码以实现基于网格搜索的初值优化（这里仅作简单示意，实际应用时需根据实际情况设定搜索范围和步长）：

% 定义搜索范围和步长
x1_range = [50, 150];
x2_range = [5, 20];
x3_range = [0.05, 0.5];
step_size = [20, 5, 0.05];

% 初始化存储最优解和对应的目标函数值变量
best_x = [];
best_residual = Inf;

for x1 = x1_range(1):step_size(1):x1_range(2)
    for x2 = x2_range(1):step_size(2):x2_range(2)
        for x3 = x3_range(1):step_size(3):x3_range(2)
            % 构造当前初值
            x0_current = [x1, x2, x3];
            
            % 使用当前初值进行拟合
            x_current = lsqcurvefit(fun, x0_current, t, y);
            
            % 计算当前解的目标函数值
            residual_current = sum((fun(x_current, t) - y).^2);
            
            % 更新最优解和目标函数值
            if residual_current < best_residual
                best_residual = residual_current;
                best_x = x_current;
            end
        end
    end
end

% 使用最优初值重新拟合
x_final = lsqcurvefit(fun, best_x, t, y);

% 绘制拟合曲线及计算评估指标
% ...（此处省略了绘图和计算R2、RMSE的部分，同原代码）

以上代码将遍历指定范围内的所有初值组合，寻找最佳拟合参数。请注意，这将显著增加计算量，特别是当搜索空间很大时，因此在实际操作时应合理设置搜索范围和步长。