Mercurial > hg > octave-max

comparison scripts/polynomial/polyfit.m @ 3091:b06dcbb6b3b1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1997-10-10 16:26:47 by jwe]
author jwe
date Fri, 10 Oct 1997 16:26:48 +0000
parents db6d57d718f7
children 4bb976b250bf
comparison
equal deleted inserted replaced
3090:63bda47c6512 3091:b06dcbb6b3b1
1 ## Copyright (C) 1996, 1997 John W. Eaton 1 ## Copyright (C) 1996 John W. Eaton
2 ## 2 ##
3 ## This file is part of Octave. 3 ## This file is part of Octave.
4 ## 4 ##
5 ## Octave is free software; you can redistribute it and/or modify it 5 ## Octave is free software; you can redistribute it and/or modify it
6 ## under the terms of the GNU General Public License as published by 6 ## under the terms of the GNU General Public License as published by
43 43
44 if (! (is_scalar (n) && n >= 0 && ! isinf (n) && n == round (n))) 44 if (! (is_scalar (n) && n >= 0 && ! isinf (n) && n == round (n)))
45 error ("polyfit: n must be a nonnegative integer"); 45 error ("polyfit: n must be a nonnegative integer");
46 endif 46 endif
47 47
48 y_is_row_vector = (rows (y) == 1);
49
48 l = length (x); 50 l = length (x);
49 x = reshape (x, l, 1); 51 x = reshape (x, l, 1);
50 y = reshape (y, l, 1); 52 y = reshape (y, l, 1);
51 53
52 ## Unfortunately, the economy QR factorization doesn't really save 54 ## Unfortunately, the economy QR factorization doesn't really save
59 ## XXX FIXME XXX -- this is probably not so good for extreme values of 61 ## XXX FIXME XXX -- this is probably not so good for extreme values of
60 ## N or X... 62 ## N or X...
61 63
62 X = (x * ones (1, n+1)) .^ (ones (l, 1) * (0 : n)); 64 X = (x * ones (1, n+1)) .^ (ones (l, 1) * (0 : n));
63 65
64 p = flipud ((X' * X) \ (X' * y)); 66 p = (X' * X) \ (X' * y);
67
68 if (nargout == 2)
69 yf = X * p;
70 endif
71
72 p = flipud (p);
65 73
66 if (! prefer_column_vectors) 74 if (! prefer_column_vectors)
67 p = p'; 75 p = p';
68 endif 76 endif
69 77
70 if (nargout == 2) 78 if (y_is_row_vector)
71 yf = X * p; 79 yf = yf';
72 endif 80 endif
73 81
74 endfunction 82 endfunction