annotate scripts/polynomial/polyfit.m @ 3784:097af73b5f4f

[project @ 2001-02-07 17:17:13 by jwe]
author jwe
date Wed, 07 Feb 2001 17:17:30 +0000
parents f8dde1807dee
children 22bd65326ec1
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
3092
4bb976b250bf [project @ 1997-10-17 07:48:48 by jwe]
jwe
parents: 3091
diff changeset
1 ## Copyright (C) 1996, 1997 John W. Eaton
2313
5ca126254d15 [project @ 1996-07-11 21:25:22 by jwe]
jwe
parents: 2312
diff changeset
2 ##
5ca126254d15 [project @ 1996-07-11 21:25:22 by jwe]
jwe
parents: 2312
diff changeset
3 ## This file is part of Octave.
5ca126254d15 [project @ 1996-07-11 21:25:22 by jwe]
jwe
parents: 2312
diff changeset
4 ##
5ca126254d15 [project @ 1996-07-11 21:25:22 by jwe]
jwe
parents: 2312
diff changeset
5 ## Octave is free software; you can redistribute it and/or modify it
5ca126254d15 [project @ 1996-07-11 21:25:22 by jwe]
jwe
parents: 2312
diff changeset
6 ## under the terms of the GNU General Public License as published by
5ca126254d15 [project @ 1996-07-11 21:25:22 by jwe]
jwe
parents: 2312
diff changeset
7 ## the Free Software Foundation; either version 2, or (at your option)
5ca126254d15 [project @ 1996-07-11 21:25:22 by jwe]
jwe
parents: 2312
diff changeset
8 ## any later version.
5ca126254d15 [project @ 1996-07-11 21:25:22 by jwe]
jwe
parents: 2312
diff changeset
9 ##
5ca126254d15 [project @ 1996-07-11 21:25:22 by jwe]
jwe
parents: 2312
diff changeset
10 ## Octave is distributed in the hope that it will be useful, but
5ca126254d15 [project @ 1996-07-11 21:25:22 by jwe]
jwe
parents: 2312
diff changeset
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
5ca126254d15 [project @ 1996-07-11 21:25:22 by jwe]
jwe
parents: 2312
diff changeset
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
5ca126254d15 [project @ 1996-07-11 21:25:22 by jwe]
jwe
parents: 2312
diff changeset
13 ## General Public License for more details.
5ca126254d15 [project @ 1996-07-11 21:25:22 by jwe]
jwe
parents: 2312
diff changeset
14 ##
5ca126254d15 [project @ 1996-07-11 21:25:22 by jwe]
jwe
parents: 2312
diff changeset
15 ## You should have received a copy of the GNU General Public License
5ca126254d15 [project @ 1996-07-11 21:25:22 by jwe]
jwe
parents: 2312
diff changeset
16 ## along with Octave; see the file COPYING. If not, write to the Free
5ca126254d15 [project @ 1996-07-11 21:25:22 by jwe]
jwe
parents: 2312
diff changeset
17 ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
5ca126254d15 [project @ 1996-07-11 21:25:22 by jwe]
jwe
parents: 2312
diff changeset
18 ## 02111-1307, USA.
2303
5cffc4b8de57 [project @ 1996-06-24 09:15:24 by jwe]
jwe
parents: 2261
diff changeset
19
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3246
diff changeset
20 ## -*- texinfo -*-
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3246
diff changeset
21 ## @deftypefn {Function File} {[@var{p}, @var{yf}] =} polyfit (@var{x}, @var{y}, @var{n})
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3246
diff changeset
22 ## Return the coefficients of a polynomial @var{p}(@var{x}) of degree
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3418
diff changeset
23 ## @var{n} that minimizes
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3246
diff changeset
24 ## @iftex
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3246
diff changeset
25 ## @tex
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3246
diff changeset
26 ## $$
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3246
diff changeset
27 ## \sum_{i=1}^N (p(x_i) - y_i)^2
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3246
diff changeset
28 ## $$
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3246
diff changeset
29 ## @end tex
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3246
diff changeset
30 ## @end iftex
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3246
diff changeset
31 ## @ifinfo
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3246
diff changeset
32 ## @code{sumsq (p(x(i)) - y(i))},
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3246
diff changeset
33 ## @end ifinfo
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3246
diff changeset
34 ## to best fit the data in the least squares sense.
3418
ca92c9d3f882 [project @ 2000-01-12 03:07:47 by jwe]
jwe
parents: 3368
diff changeset
35 ##
ca92c9d3f882 [project @ 2000-01-12 03:07:47 by jwe]
jwe
parents: 3368
diff changeset
36 ## The polynomial coefficients are returned in a row vector if @var{x}
ca92c9d3f882 [project @ 2000-01-12 03:07:47 by jwe]
jwe
parents: 3368
diff changeset
37 ## and @var{y} are both row vectors; otherwise, they are returned in a
ca92c9d3f882 [project @ 2000-01-12 03:07:47 by jwe]
jwe
parents: 3368
diff changeset
38 ## column vector.
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3418
diff changeset
39 ##
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3246
diff changeset
40 ## If two output arguments are requested, the second contains the values of
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3246
diff changeset
41 ## the polynomial for each value of @var{x}.
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3246
diff changeset
42 ## @end deftypefn
2311
2b5788792cad [project @ 1996-07-11 20:18:38 by jwe]
jwe
parents: 2303
diff changeset
43
2312
204cc7db6f4a [project @ 1996-07-11 21:20:36 by jwe]
jwe
parents: 2311
diff changeset
44 ## Author: KH <Kurt.Hornik@ci.tuwien.ac.at>
204cc7db6f4a [project @ 1996-07-11 21:20:36 by jwe]
jwe
parents: 2311
diff changeset
45 ## Created: 13 December 1994
204cc7db6f4a [project @ 1996-07-11 21:20:36 by jwe]
jwe
parents: 2311
diff changeset
46 ## Adapted-By: jwe
204cc7db6f4a [project @ 1996-07-11 21:20:36 by jwe]
jwe
parents: 2311
diff changeset
47
3061
1a5fe3010f09 [project @ 1997-06-25 18:36:19 by jwe]
jwe
parents: 2847
diff changeset
48 function [p, yf] = polyfit (x, y, n)
2325
b5568c31ee2c [project @ 1996-07-15 22:20:21 by jwe]
jwe
parents: 2313
diff changeset
49
b5568c31ee2c [project @ 1996-07-15 22:20:21 by jwe]
jwe
parents: 2313
diff changeset
50
2261
1b6e1629fb91 [project @ 1996-05-23 00:52:07 by jwe]
jwe
parents:
diff changeset
51 if (nargin != 3)
1b6e1629fb91 [project @ 1996-05-23 00:52:07 by jwe]
jwe
parents:
diff changeset
52 usage ("polyfit (x, y, n)");
1b6e1629fb91 [project @ 1996-05-23 00:52:07 by jwe]
jwe
parents:
diff changeset
53 endif
2325
b5568c31ee2c [project @ 1996-07-15 22:20:21 by jwe]
jwe
parents: 2313
diff changeset
54
2261
1b6e1629fb91 [project @ 1996-05-23 00:52:07 by jwe]
jwe
parents:
diff changeset
55 if (! (is_vector (x) && is_vector (y) && size (x) == size (y)))
1b6e1629fb91 [project @ 1996-05-23 00:52:07 by jwe]
jwe
parents:
diff changeset
56 error ("polyfit: x and y must be vectors of the same size");
1b6e1629fb91 [project @ 1996-05-23 00:52:07 by jwe]
jwe
parents:
diff changeset
57 endif
2325
b5568c31ee2c [project @ 1996-07-15 22:20:21 by jwe]
jwe
parents: 2313
diff changeset
58
2261
1b6e1629fb91 [project @ 1996-05-23 00:52:07 by jwe]
jwe
parents:
diff changeset
59 if (! (is_scalar (n) && n >= 0 && ! isinf (n) && n == round (n)))
1b6e1629fb91 [project @ 1996-05-23 00:52:07 by jwe]
jwe
parents:
diff changeset
60 error ("polyfit: n must be a nonnegative integer");
1b6e1629fb91 [project @ 1996-05-23 00:52:07 by jwe]
jwe
parents:
diff changeset
61 endif
2325
b5568c31ee2c [project @ 1996-07-15 22:20:21 by jwe]
jwe
parents: 2313
diff changeset
62
3091
b06dcbb6b3b1 [project @ 1997-10-10 16:26:47 by jwe]
jwe
parents: 3069
diff changeset
63 y_is_row_vector = (rows (y) == 1);
b06dcbb6b3b1 [project @ 1997-10-10 16:26:47 by jwe]
jwe
parents: 3069
diff changeset
64
2261
1b6e1629fb91 [project @ 1996-05-23 00:52:07 by jwe]
jwe
parents:
diff changeset
65 l = length (x);
1b6e1629fb91 [project @ 1996-05-23 00:52:07 by jwe]
jwe
parents:
diff changeset
66 x = reshape (x, l, 1);
1b6e1629fb91 [project @ 1996-05-23 00:52:07 by jwe]
jwe
parents:
diff changeset
67 y = reshape (y, l, 1);
2325
b5568c31ee2c [project @ 1996-07-15 22:20:21 by jwe]
jwe
parents: 2313
diff changeset
68
3061
1a5fe3010f09 [project @ 1997-06-25 18:36:19 by jwe]
jwe
parents: 2847
diff changeset
69 X = (x * ones (1, n+1)) .^ (ones (l, 1) * (0 : n));
2261
1b6e1629fb91 [project @ 1996-05-23 00:52:07 by jwe]
jwe
parents:
diff changeset
70
3190
70eb3f4136cd [project @ 1998-10-19 22:37:56 by jwe]
jwe
parents: 3105
diff changeset
71 p = X \ y;
3091
b06dcbb6b3b1 [project @ 1997-10-10 16:26:47 by jwe]
jwe
parents: 3069
diff changeset
72
b06dcbb6b3b1 [project @ 1997-10-10 16:26:47 by jwe]
jwe
parents: 3069
diff changeset
73 if (nargout == 2)
b06dcbb6b3b1 [project @ 1997-10-10 16:26:47 by jwe]
jwe
parents: 3069
diff changeset
74 yf = X * p;
3105
f936c7f5074f [project @ 1997-11-19 04:36:17 by jwe]
jwe
parents: 3092
diff changeset
75
f936c7f5074f [project @ 1997-11-19 04:36:17 by jwe]
jwe
parents: 3092
diff changeset
76 if (y_is_row_vector)
3246
a41cc560087a [project @ 1999-06-19 06:46:20 by jwe]
jwe
parents: 3238
diff changeset
77 yf = yf.';
3105
f936c7f5074f [project @ 1997-11-19 04:36:17 by jwe]
jwe
parents: 3092
diff changeset
78 endif
3091
b06dcbb6b3b1 [project @ 1997-10-10 16:26:47 by jwe]
jwe
parents: 3069
diff changeset
79 endif
b06dcbb6b3b1 [project @ 1997-10-10 16:26:47 by jwe]
jwe
parents: 3069
diff changeset
80
b06dcbb6b3b1 [project @ 1997-10-10 16:26:47 by jwe]
jwe
parents: 3069
diff changeset
81 p = flipud (p);
2261
1b6e1629fb91 [project @ 1996-05-23 00:52:07 by jwe]
jwe
parents:
diff changeset
82
3418
ca92c9d3f882 [project @ 2000-01-12 03:07:47 by jwe]
jwe
parents: 3368
diff changeset
83 if (y_is_row_vector && rows (x) == 1)
ca92c9d3f882 [project @ 2000-01-12 03:07:47 by jwe]
jwe
parents: 3368
diff changeset
84 p = p';
2261
1b6e1629fb91 [project @ 1996-05-23 00:52:07 by jwe]
jwe
parents:
diff changeset
85 endif
1b6e1629fb91 [project @ 1996-05-23 00:52:07 by jwe]
jwe
parents:
diff changeset
86
1b6e1629fb91 [project @ 1996-05-23 00:52:07 by jwe]
jwe
parents:
diff changeset
87 endfunction