1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
|
function [y, delta] = polyval(p,x,S,mu)
//y = polyval(p,x) returns the value of a polynomial of degree n evaluated at x. The input argument p is a vector of length n+1 whose elements are the coefficients in descending powers of the polynomial to be evaluated.
//x --can be a matrix or a vector. In either case, polyval evaluates p at each element of x.
//[y,delta] = polyval(p,x,S) uses the optional output structure S generated by polyfit to generate error estimates delta.
//delta --is an estimate of the standard deviation of the error in predicting a future observation at x by p(x).
// If the coefficients in p are least squares estimates computed by polyfit, and the errors in the data input to poly// if it are independent, normal, and have constant variance, then y±delta contains at least 50% of the
// predictions of future observations at x.
//y = polyval(p,x,[],mu) or [y,delta] = polyval(p,x,S,mu) use ˆx=(x−μ1)/μ2 in place of x. In this equation, μ1=mean(x) // and μ2=std(x). The centering and scaling parameters mu = [μ1,μ2] is optional
//EXAMPLES:
//p = [3 2 1];
//y=polyval(p,[5 7 9])
//EXPECTED OUTPUT:
//y= 86 162 262
if ~(isvector(p) | isempty(p)) // Check input is a vector
error('polyval:InvalidP');
end
nc = length(p);
if isscalar(x) & (argn(2) < 3) & nc>0 & (abs(x)<%inf) & and(abs(p(:))<%inf)
// Make it scream for scalar x. Polynomial evaluation can be
// implemented as a recursive digital filter.
y = filter(1,[1 -real(x)],p);
y = y(nc);
return
end
siz_x = size(x);
if argn(2) == 4
x = (x - mu(1))/mu(2);
end
// Use Horner's method for general case where X is an array.
y = zeros(size(x,1),size(x,2));
if nc>0, y(:) = p(1); end
for i=2:nc
y = x .* y + p(i);
end
if argn(1) > 1
if argn(2) < 3 | isempty(S)
error('polyval:Requires S');
end
// Extract parameters from S
if isstruct(S), // Use output structure from polyfit.
R = S.R;
df = S.df;
normr = S.normr;
else // Use output matrix from previous versions of polyfit.
[ms,ns] = size(S);
if (ms ~= ns+2) | (nc ~= ns)
error('polyval:SizeS');
end
R = S(1:nc,1:nc);
df = S(nc+1,1);
normr = S(nc+2,1);
end
// Construct Vandermonde matrix for the new X.
x = x(:);
V(:,nc) = ones(length(x),1,class(x));
for j = nc-1:-1:1
V(:,j) = x.*V(:,j+1);
end
// S is a structure containing three elements: the triangular factor of
// the Vandermonde matrix for the original X, the degrees of freedom,
// and the norm of the residuals.
E = V/R;
e = sqrt(1+sum(E.*E,2));
if df == 0
warning('polyval:ZeroDOF');
delta = %inf(size(e));
else
delta = normr/sqrt(df)*e;
end
delta = matrix(delta,siz_x);
end
endfunction
|
