Mercurial > hg > octave-lyh

comparison doc/interpreter/poly.txi @ 14509:a88f8e4fae56

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
New Function, splinefit.m * __splinefit__.m: New private file. Jonas Lundgren's splinefit.m with BSD License. Jonas emailed this version to Octave's developers. * splinefit.m: New File. Piece-wise polynomial fit. This is a wrapper for __splinefit__.m. The wrapper allows for Octave's tex-info documentation, demos, and tests to be added. In addition the input syntax has been sligtly modified to allow new options to be added without breaking compatiblity. * doc/splineimages.m: New file to produce splineimages<#> for the docs. * doc/images: Include splineimages.m and the figues for the docs. * scripts/polynomial/module.mk: Add new files. * doc/interpreter/poly.txi: Minimal description of splinefit.
author Ben Abbott <bpabbott@mac.com>
date Thu, 29 Mar 2012 19:13:21 -0400
parents 72c96de7a403
children 8985cbbd2fe4
comparison
equal deleted inserted replaced
14508:0901f926ed50 14509:a88f8e4fae56
130 @section Polynomial Interpolation 130 @section Polynomial Interpolation
131 131
132 Octave comes with good support for various kinds of interpolation, 132 Octave comes with good support for various kinds of interpolation,
133 most of which are described in @ref{Interpolation}. One simple alternative 133 most of which are described in @ref{Interpolation}. One simple alternative
134 to the functions described in the aforementioned chapter, is to fit 134 to the functions described in the aforementioned chapter, is to fit
135 a single polynomial to some given data points. To avoid a highly 135 a single polynomial, or a piecewise polynomial (spline) to some given
136 fluctuating polynomial, one most often wants to fit a low-order polynomial 136 data points. To avoid a highly fluctuating polynomial, one most often
137 to data. This usually means that it is necessary to fit the polynomial 137 wants to fit a low-order polynomial to data. This usually means that it
138 in a least-squares sense, which just is what the @code{polyfit} function does. 138 is necessary to fit the polynomial in a least-squares sense, which just
139 is what the @code{polyfit} function does.
139 140
140 @DOCSTRING(polyfit) 141 @DOCSTRING(polyfit)
141 142
142 In situations where a single polynomial isn't good enough, a solution 143 In situations where a single polynomial isn't good enough, a solution
143 is to use several polynomials pieced together. The function @code{mkpp} 144 is to use several polynomials pieced together. The function
144 creates a piecewise polynomial, @code{ppval} evaluates the function 145 @code{splinefit} fits a peicewise polynomial (spline) to a set of
145 created by @code{mkpp}, and @code{unmkpp} returns detailed information 146 data.
146 about the function. 147
148 @DOCSTRING(splinefit)
149
150 The number of @var{breaks} (or knots) used to construct the piecewise
151 polynomial is a significant factor in suppressing the noise present in
152 the input data, @var{x} and @var{y}. This is demostrated by the example
153 below.
154
155 @example
156 @group
157 x = 2 * pi * rand (1, 200);
158 y = sin (x) + sin (2 * x) + 0.2 * randn (size (x));
159 ## Uniform breaks
160 breaks = linspace (0, 2 * pi, 41); % 41 breaks, 40 pieces
161 pp1 = splinefit (x, y, breaks);
162 ## Breaks interpolated from data
163 pp2 = splinefit (x, y, 10); % 11 breaks, 10 pieces
164 ## Plot
165 xx = linspace (0, 2 * pi, 400);
166 y1 = ppval (pp1, xx);
167 y2 = ppval (pp2, xx);
168 plot (x, y, ".", xx, [y1; y2])
169 axis tight
170 ylim auto
171 legend (@{"data", "41 breaks, 40 pieces", "11 breaks, 10 pieces"@})
172 @end group
173 @end example
174
175 @ifnotinfo
176 @noindent
177 The result of which can be seen in @ref{fig:splinefit1}.
178
179 @float Figure,fig:splinefit1
180 @center @image{splinefit1,4in}
181 @caption{Comparison of a fitting a piecewise polynomial with 41 breaks to one
182 with 11 breaks. The fit with the large number of breaks exhibits a fast ripple
183 that is not present in the underlying function.}
184 @end float
185 @end ifnotinfo
186
187 The piece-wise polynomial fit provided by @code{splinefit} provides
188 continuous derivatives up to the @var{order} of the fit. This can
189 be demonstrated by the code
190
191 @example
192 @group
193 ## Data (200 points)
194 x = 2 * pi * rand (1, 200);
195 y = sin (x) + sin (2 * x) + 0.1 * randn (size (x));
196 ## Piecewise constant
197 pp1 = splinefit (x, y, 8, "order", 0);
198 ## Piecewise linear
199 pp2 = splinefit (x, y, 8, "order", 1);
200 ## Piecewise quadratic
201 pp3 = splinefit (x, y, 8, "order", 2);
202 ## Piecewise cubic
203 pp4 = splinefit (x, y, 8, "order", 3);
204 ## Piecewise quartic
205 pp5 = splinefit (x, y, 8, "order", 4);
206 ## Plot
207 xx = linspace (0, 2 * pi, 400);
208 y1 = ppval (pp1, xx);
209 y2 = ppval (pp2, xx);
210 y3 = ppval (pp3, xx);
211 y4 = ppval (pp4, xx);
212 y5 = ppval (pp5, xx);
213 plot (x, y, ".", xx, [y1; y2; y3; y4; y5])
214 axis tight
215 ylim auto
216 legend (@{"data", "order 0", "order 1", "order 2", "order 3", "order 4"@})
217 @end group
218 @end example
219
220 @ifnotinfo
221 @noindent
222 The result of which can be seen in @ref{fig:splinefit2}.
223
224 @float Figure,fig:splinefit2
225 @center @image{splinefit2,4in}
226 @caption{Comparison of a piecewise constant, linear, quadratic, cubic, and
227 quartic polynomials with 8 breaks to noisey data. The higher order solutions
228 more accurately represent the underlying function, but come with the
229 expense of computational complexity.}
230 @end float
231 @end ifnotinfo
232
233 When the underlying function to provide a fit to is periodice, @code{splitfit}
234 is able to apply the boundary conditions needed to manifest a periodic fit.
235 This is demonstrated by the code below.
236
237 @example
238 @group
239 ## Data (100 points)
240 x = 2 * pi * [0, (rand (1, 98)), 1];
241 y = sin (x) - cos (2 * x) + 0.2 * randn (size (x));
242 ## No constraints
243 pp1 = splinefit (x, y, 10, "order", 5);
244 ## Periodic boundaries
245 pp2 = splinefit (x, y, 10, "order", 5, "periodic", true);
246 ## Plot
247 xx = linspace (0, 2 * pi, 400);
248 y1 = ppval (pp1, xx);
249 y2 = ppval (pp2, xx);
250 plot (x, y, ".", xx, [y1; y2])
251 axis tight
252 ylim auto
253 legend (@{"data", "no constraints", "periodic"@})
254 @end group
255 @end example
256
257 @ifnotinfo
258 @noindent
259 The result of which can be seen in @ref{fig:splinefit3}.
260
261 @float Figure,fig:splinefit3
262 @center @image{splinefit3,4in}
263 @caption{Comparison of piecewise cubic fits to a noisy periodic function with,
264 and without, periodic boundary conditions.}
265 @end float
266 @end ifnotinfo
267
268 More complex constaints may be added as well. For example, the code below
269 illustrates a periodic fit with values that have been clamped end point values
270 and a second periodic fit wihh hinged end point values.
271
272 @example
273 @group
274 ## Data (200 points)
275 x = 2 * pi * rand (1, 200);
276 y = sin (2 * x) + 0.1 * randn (size (x));
277 ## Breaks
278 breaks = linspace (0, 2 * pi, 10);
279 ## Clamped endpoints, y = y" = 0
280 xc = [0, 0, 2*pi, 2*pi];
281 cc = [(eye (2)), (eye (2))];
282 con = struct ("xc", xc, "cc", cc);
283 pp1 = splinefit (x, y, breaks, "constraints", con);
284 ## Hinged periodic endpoints, y = 0
285 con = struct ("xc", 0);
286 pp2 = splinefit (x, y, breaks, "constraints", con, "periodic", true);
287 ## Plot
288 xx = linspace (0, 2 * pi, 400);
289 y1 = ppval (pp1, xx);
290 y2 = ppval (pp2, xx);
291 plot (x, y, ".", xx, [y1; y2])
292 axis tight
293 ylim auto
294 legend (@{"data", "clamped", "hinged periodic"@})
295 @end group
296 @end example
297
298 @ifnotinfo
299 @noindent
300 The result of which can be seen in @ref{fig:splinefit4}.
301
302 @float Figure,fig:splinefit4
303 @center @image{splinefit4,4in}
304 @caption{Comparison of two periodic piecewise cubic fits to a noisy periodic
305 signal. One fit has its end points clamped and the second has its end points
306 hinged.}
307 @end float
308 @end ifnotinfo
309
310 The @code{splinefit} function also provides the convenience of a @var{robust}
311 fitting, where the effect of outlying data is reduced. In the example below,
312 three different fits are provided. Two with differing levels of outlier
313 suppressin and a third illustrating the non-robust solution.
314
315 @example
316 @group
317 ## Data
318 x = linspace (0, 2*pi, 200);
319 y = sin (x) + sin (2 * x) + 0.05 * randn (size (x));
320 ## Add outliers
321 x = [x, linspace(0,2*pi,60)];
322 y = [y, -ones(1,60)];
323 ## Fit splines with hinged conditions
324 con = struct ("xc", [0, 2*pi]);
325 ## Robust fitting, beta = 0.25
326 pp1 = splinefit (x, y, 8, "constraints", con, "beta", 0.25);
327 ## Robust fitting, beta = 0.75
328 pp2 = splinefit (x, y, 8, "constraints", con, "beta", 0.75);
329 ## No robust fitting
330 pp3 = splinefit (x, y, 8, "constraints", con);
331 ## Plot
332 xx = linspace (0, 2*pi, 400);
333 y1 = ppval (pp1, xx);
334 y2 = ppval (pp2, xx);
335 y3 = ppval (pp3, xx);
336 plot (x, y, ".", xx, [y1; y2; y3])
337 legend (@{"data with outliers","robust, beta = 0.25", ...
338 "robust, beta = 0.75", "no robust fitting"@})
339 axis tight
340 ylim auto
341 @end group
342 @end example
343
344 @ifnotinfo
345 @noindent
346 The result of which can be seen in @ref{fig:splinefit6}.
347
348 @float Figure,fig:splinefit6
349 @center @image{splinefit6,4in}
350 @caption{Comparison of two differnet levels of robust fitting (@var{beta} = 0.25 and 0.75) to noisy data combined with outlying data. A standard fit is also included.}
351 @end float
352 @end ifnotinfo
353
354 The function, @code{ppval}, evaluates the piecewise polyomials, created
355 by @code{mkpp} or other means, and @code{unmkpp} returns detailed
356 information about the piecewise polynomial.
147 357
148 The following example shows how to combine two linear functions and a 358 The following example shows how to combine two linear functions and a
149 quadratic into one function. Each of these functions is expressed 359 quadratic into one function. Each of these functions is expressed
150 on adjoined intervals. 360 on adjoined intervals.
151 361