Matlab正弦拟合程序中参数误差与协方差的实现问询

我编写了一段Matlab代码实现数据的正弦拟合，代码如下：

data = importdata('analisipicco.txt') ;  
x = data(:,1) ; y = data(:,2) ; 
yu = max(y);
yl = min(y);
yr = (yu-yl);                               % Range of ‘y’
yz = y-yu+(yr/2);
zx = x(yz .* circshift(yz,1) <= 0);     % Find zero-crossings
per = 2*mean(diff(zx));                     % Estimate period
ym = mean(y);                               % Estimate offset

fit = @(b,x)  b(1).*(sin(2*pi*x./b(2) + 2*pi/b(3))) + b(4);    % Function to fit
fcn = @(b) sum((fit(b,x) - y).^2);                              % Least-Squares cost function
s = fminsearch(fcn, [yr;  per;  -1;  ym]);                      % Minimise Least-Squares

xp = linspace(min(x),max(x));

xlabel("passi motore");
ylabel("intensità (u.a.)");
figure(1)
plot(x,y,'o',  xp,fit(s,xp), 'r')
grid

输出参数向量s（对应函数中的b）各元素含义为：

• s(1)：正弦波振幅（y单位）
• s(2)：周期（x单位）
• s(3)：相位（相位值为s(2)/(2*s(3))，x单位）
• s(4)：偏移量（y单位）

我的数据集如下：

-200 183966
-192 189734
-184 195724
-176 201663
-168 207557
-160 213278
-152 219000
-144 224677
-136 229500
-128 236024
-120 241968
-112 247787
-104 252963
-96 257491
-88 261967
-80 267373
-72 273494
-64 278599
-56 281476
-48 282610
-40 283097
-32 283839
-24 284971
-16 286169
-8  287164
0   287968
8   288561
16  288626
24  288107
32  286967
40  285132
48  282828
56  279847
64  276296
72  272299
80  268080
88  263564
96  258926
104 254052
112 248894
120 243694
128 238177
136 232665
144 227143
152 221959
160 216874
168 211678
176 206540
184 201537
192 196748
200 192091

程序运行正常，但我想求解拟合参数的误差及协方差矩阵。我知道可以用矩阵法实现，但作为Matlab新手，不知道如何编写代码，想咨询是否有内置函数可直接计算参数误差？

Matlab里确实有便捷的内置工具可以计算拟合参数的误差和协方差矩阵，推荐两种常用方法：

方法一：使用lsqcurvefit（优化工具箱） #

该函数完成最小二乘拟合的同时，可输出雅克比矩阵，进而推导协方差矩阵和参数误差。修改后的代码如下：

data = importdata('analisipicco.txt') ;  
x = data(:,1) ; y = data(:,2) ; 
yu = max(y);
yl = min(y);
yr = (yu-yl);                               % Range of ‘y’
yz = y-yu+(yr/2);
zx = x(yz .* circshift(yz,1) <= 0);     % Find zero-crossings
per = 2*mean(diff(zx));                     % Estimate period
ym = mean(y);                               % Estimate offset

% 定义拟合函数
fit_fun = @(b,x) b(1).*sin(2*pi*x./b(2) + 2*pi/b(3)) + b(4);

% 初始参数猜测
initial_guess = [yr; per; -1; ym];

% 调用lsqcurvefit，获取雅克比矩阵
options = optimoptions('lsqcurvefit','Display','off');
[s, resnorm, ~, ~, ~, ~, jacobian] = lsqcurvefit(fit_fun, initial_guess, x, y, [], [], options);

% 计算协方差矩阵与参数标准误差
n = length(y);
p = length(s);
sigma_sq = resnorm/(n-p);  % 残差方差估计
cov_matrix = sigma_sq * inv(jacobian'*jacobian);
param_errors = sqrt(diag(cov_matrix));  % 参数标准误差

% 输出结果
disp('拟合参数：');
disp(s);
disp('参数标准误差：');
disp(param_errors);
disp('协方差矩阵：');
disp(cov_matrix);

% 绘图
xp = linspace(min(x),max(x));
xlabel("passi motore");
ylabel("intensità (u.a.)");
figure(1)
plot(x,y,'o', xp, fit_fun(s,xp), 'r')
grid

方法二：使用曲线拟合工具箱的fit函数 #

若安装了曲线拟合工具箱，代码更直观，可直接获取置信区间和协方差矩阵：

data = importdata('analisipicco.txt') ;  
x = data(:,1) ; y = data(:,2) ; 
yu = max(y);
yl = min(y);
yr = (yu-yl);                               % Range of ‘y’
yz = y-yu+(yr/2);
zx = x(yz .* circshift(yz,1) <= 0);     % Find zero-crossings
per = 2*mean(diff(zx));                     % Estimate period
ym = mean(y);                               % Estimate offset

% 自定义拟合类型
ft = fittype('b1*sin(2*pi*x/b2 + 2*pi/b3) + b4', 'independent', 'x', 'coefficients', {'b1','b2','b3','b4'});

% 设置初始参数
opts = fitoptions(ft);
opts.StartPoint = [yr, per, -1, ym];

% 执行拟合
[fitresult, ~] = fit(x, y, ft, opts);

% 输出参数及置信区间
disp(fitresult);
% 提取协方差矩阵
cov_matrix = coeffcov(fitresult);
disp('协方差矩阵：');
disp(cov_matrix);

% 绘图
xp = linspace(min(x),max(x));
xlabel("passi motore");
ylabel("intensità (u.a.)");
figure(1)
plot(x,y,'o', xp, fitresult(xp), 'r')
grid

关键说明 #

1. lsqcurvefit通过雅克比矩阵计算协方差矩阵，sigma_sq是基于自由度的残差方差估计，参数标准误差为协方差矩阵对角线元素的平方根。
2. 曲线拟合工具箱的fit函数会自动展示参数的95%置信区间，coeffcov函数可直接提取协方差矩阵。

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