如何从多项式拟合中绘制参数误差 - MATLAB

Question

我有一个数据列表,我试图适应多项式,我试图绘制参数的95%置信区间(在Matlab中).如果我的数据是x和y

f=fit(x,y,'poly2')
plot(f,x,y)
ci=confint(f,0.95);
a_ci=ci(1,:);
b_ci=ci(2,:);

我不知道如何在此之后继续获取我的数据的最小和最大波段.有谁知道这是怎么做到的吗？

Answer 1

我可以看到你安装了曲线拟合工具箱,这很好,因为你需要它才能使下面的代码工作.

示例数据的基本拟合

让我们定义一些示例数据和一个可能的拟合函数.(我也可以poly2在这里使用,但我想保持它更一般.)

xdata = (0:0.1:1)';               % column vector!
noise = 0.1*randn(size(xdata));
ydata = xdata.^2 + noise;
f = fittype('a*x.^2 + b'); 
fit1 = fit(xdata, ydata, f, 'StartPoint', [1,1])
plot(fit1, xdata, ydata)

旁注:plot()不是我们通常的绘图功能,而是cfit-object fit1的一种方法.

拟合参数的置信区间

我们的配合使用数据来确定系数a,b就底层模型f(x)=ax2+b.你已经这样做了,但为了完整起见,你可以在这里读出任何置信区间系数的不确定性.系数按字母顺序排列,这就是我可以使用ci(1,:)的原因a,依此类推.

names = coeffnames(fit1)   % check the coefficient order!
ci = confint(fit1, 0.95);  % 2 sigma interval
a_ci = ci(1,:)
b_ci = ci(2,:)

默认情况下,Matlab使用2σ(0.95)置信区间.有些人(物理学家)更喜欢引用1σ(0.68)区间.

信心和预测带

在数据周围绘制置信带或预测带是一个好习惯- 特别是当系数相关时!但是你应该花一点时间思考你想要绘制的两个中的哪一个:

• 预测范围:如果我采用新的测量值,我认为它应该在哪里？在Matlab术语中,这被称为"观察带".
• 信心带:我认为真正的价值在哪里？在Matlab术语中,这被称为"功能带".

与系数的置信区间一样,Matlab默认使用2σ波段,我们中的物理学家将其切换为1σ区间.就其本质而言,预测带更大,因为它是模型误差(置信带!)和测量误差的组合.

还有另一个目标,一个我不完全理解的目标.Matlab和维基百科都做出了这种区分.

• Pointwise:单个测量/真值的预测/置信带有多大？在我能想到的几乎所有情况下,这都是你想要作为一名物理学家所要求的.
• 同时:如果您想要一组具有给定置信度的所有新测量值/所有预测点位于波段内,您有多大的预测/置信带？

在我个人看来,"同步乐队"不是乐队!对于具有n个点的测量,它应该是n个单独的误差棒!

预测/置信区分和逐点/同时区分为图表周围的"带"提供了总共四个选项.Matlab使2σ 逐点预测带易于访问,但您似乎感兴趣的是2σ 逐点置信带.这是一个比较繁琐一点阴谋,因为你必须指定在其评估预测带伪数据:

x_dummy = linspace(min(xdata), max(xdata), 100);
figure(1); clf(1);
hold all
plot(xdata,ydata,'.')
plot(fit1)    % by default, evaluates the fit over the currnet XLim
% use "functional" (confidence!) band; use "simultaneous"=off
conf1 = predint(fit1,x_dummy,0.95,'functional','off');
plot(x_dummy, conf1, 'r--')
hold off

请注意,置信带x=0等于拟合系数的置信区间b!

外推

如果您想要推断到数据范围未涵盖的x值,您可以评估拟合和预测/置信区间的更大范围:

x_range = [0, 2];
x_dummy = linspace(x_range(1), x_range(2), 100);
figure(1); clf(1);
hold all
plot(xdata,ydata,'.')
xlim(x_range)
plot(fit1)
conf1 = predint(fit1,x_dummy,0.68,'functional','off');
plot(x_dummy, conf1, 'r--')
hold off