matlab中polyconf,Polynomial confidence intervals

灰太狼 2023-01-15 03:53 13阅读 0赞

Fit a polynomial to a sample data set, and estimate the 95% prediction intervals and the roots of the fitted polynomial. Plot the data and the estimations, and display the fitted polynomial expression using the helper function polystr, whose code appears at the end of this example.

Generate sample data points (x,y) with a quadratic trend.

rng('default') % For reproducibility

x = -5:5;

y = x.^2 - 20\*x - 3 + 5\*randn(size(x));

Fit a second degree polynomial to the data by using polyfit.

degree = 2; % Degree of the fit

\[p,S\] = polyfit(x,y,degree);

Estimate the 95% prediction intervals by using polyconf.

alpha = 0.05; % Significance level

\[yfit,delta\] = polyconf(p,x,S,'alpha',alpha);

Plot the data, fitted polynomial, and prediction intervals. Display the fitted polynomial expression using the helper function polystr.

plot(x,y,'b+')

hold on

plot(x,yfit,'g-')

plot(x,yfit-delta,'r--',x,yfit+delta,'r--')

legend('Data','Fit','95% Prediction Intervals')

title(\['Fit: ',texlabel(polystr(round(p,2)))\])

hold off

![408cfacff1b3b8191a9e8cc88829c898.png][]

Find the roots of the polynomial p.

r = roots(p)

r = 2×1

17.5152

\-0.1017

Because the roots are real values, you can plot them as well. Estimate the fitted values and prediction intervals for the x interval that includes the roots. Then, plot the roots and the estimations.

if isreal(r)

xmin = min(\[r(:);x(:)\]);

xrange = range(\[r(:);x(:)\]);

xExtended = linspace(xmin - 0.1\*xrange, xmin + 1.1\*xrange,1000);

\[yfitExtended,deltaExtended\] = polyconf(p,xExtended,S,'alpha',alpha);

plot(x,y,'b+')

hold on

plot(xExtended,yfitExtended,'g-')

plot(r,zeros(size(r)),'ko')

plot(xExtended,yfitExtended-deltaExtended,'r--')

plot(xExtended,yfitExtended+deltaExtended,'r--')

plot(xExtended,zeros(size(xExtended)),'k-')

legend('Data','Fit','Roots of Fit','95% Prediction Intervals')

title(\['Fit: ',texlabel(polystr(round(p,2)))\])

axis tight

hold off

end

![cefdbc9669bcb92cd2971dec5d82df3a.png][]

Alternatively, you can use polytool for interactive polynomial fitting.

polytool(x,y,degree,alpha)

![6b367a7c1b53e6844fd1550d357dc1ce.png][]

Helper Function

The polystr.m file defines the polystr helper function.

type polystr.m % Display contents of polystr.m file

function s = polystr(p)

% POLYSTR Converts a vector of polynomial coefficients to a character string.

% S is the string representation of P.

if all(p == 0) % All coefficients are 0.

s = '0';

else

d = length(p) - 1; % Degree of polynomial.

s = \[\]; % Initialize s.

for a = p

if a ~= 0 % Coefficient is nonzero.

if ~isempty(s) % String is not empty.

if a > 0

s = \[s ' + '\]; % Add next term.

else

s = \[s ' - '\]; % Subtract next term.

a = -a; % Value to subtract.

end

if a ~= 1 || d == 0 % Add coefficient if it is ~=1 or polynomial is constant.

s = \[s num2str(a)\];

if d > 0 % For nonconstant polynomials, add \*.

s = \[s '\*'\];

end

if d >= 2 % For terms of degree > 1, add power of x.

s = \[s 'x^' int2str(d)\];

elseif d == 1 % No power on x term.

s = \[s 'x'\];

end

d = d - 1; % Increment loop: Add term of next lowest degree.

end

[408cfacff1b3b8191a9e8cc88829c898.png]: /images/20221022/552560de887944ad92990537cfedd82e.png
[cefdbc9669bcb92cd2971dec5d82df3a.png]: /images/20221022/440ce4cd95054b8d9dfa9bc2585797a6.png
[6b367a7c1b53e6844fd1550d357dc1ce.png]: /images/20221022/3bc1ecd2692e49e78902369545dbdc56.png