1 /* 2 * $Id: fitlsq.i,v 1.1 2005-09-18 22:05:56 dhmunro Exp $ 3 * Least squares fit a piecewise linear function to data. 4 */ 5 /* Copyright (c) 2005, The Regents of the University of California. 6 * All rights reserved. 7 * This file is part of yorick (http://yorick.sourceforge.net). 8 * Read the accompanying LICENSE file for details. 9 */ 10 11 func fitlsq(y, x, xp, weight=, stats=) 12 /* DOCUMENT yp= fitlsq(y, x, xp) 13 ... 14 yfit= interp(yp, xp, xfit) 15 performs a least squares fit to the data points (X, Y). The input 16 XP are the abcissas of the piecewise linear function passing through 17 (XP, YP) which is the best fit to the data (X,Y) in a least squares 18 sense. The XP must be strictly monotone, either increasing or 19 decreasing. As for interp, the piecewise linear fit is taken to be 20 constant outside the limits of the XP. 21 22 The result YP is linearly interpolated through any consecutive 23 intervals of XP which contain no data points X, and linearly 24 extrapolated beyond the extreme values of X (if any of the intervals 25 of XP lie outside these extremes). 26 27 A WEIGHT keyword of the same length as X and Y may be supplied in 28 order to weight the various data points differently; a typical 29 WEIGHT function is 1/sigma^2 where sigma are the standard deviations 30 associated with the Y values. 31 32 SEE ALSO: interp 33 */ 34 { 35 np1= numberof(xp)+1; 36 37 /* bin the input data into the xp, and create an extended version 38 of xp for which the bins may be directly used as an index list */ 39 l= digitize(x, xp); 40 dx= (xp(0)-xp(1))*1.e30; 41 xx= grow(xp(1)-dx, xp, xp(0)+dx); 42 43 xl= xx(l); 44 xu= xx(l+1); 45 dx= xu-xl; 46 g= (x-xl)/dx; 47 h= (xu-x)/dx; 48 49 if (is_void(weight)) weight= 1.0; 50 hy= histogram(l, h*y*weight, top=np1); 51 hh= histogram(l, h*h*weight, top=np1); 52 gh= histogram(l, g*h*weight, top=np1); 53 gg= histogram(l, g*g*weight, top=np1); 54 gy= histogram(l, g*y*weight, top=np1); 55 56 diag= hh(2:0)+gg(1:-1); 57 rhs= hy(2:0)+gy(1:-1); 58 59 /* the triadiagonal system will be singular if there are any pairs 60 of consecutive bins -- remove these first */ 61 list= where(diag); 62 diag= diag(list); 63 rhs= rhs(list); 64 off= gh(list)(2:0); /* sub and super diagonal */ 65 xpp= double(xp(list)); 66 yp= TDsolve(off, diag, off, rhs); 67 68 /* special treatment if endpoints removed allows linear extrapolation 69 of the fit until the true endpoints of the xp */ 70 if (list(1)>1) { 71 yp= grow(yp(1)+(xp(1)-xpp(1))*(yp(2)-yp(1))/(xpp(2)-xpp(1)), yp); 72 xpp= grow(xp(1), xpp); 73 } 74 if (list(0)<numberof(xp)) { 75 grow, yp, yp(0)+(xp(0)-xpp(0))*(yp(0)-yp(-1))/(xpp(0)-xpp(-1)); 76 grow, xpp, xp(0) 77 } 78 79 /* add back any removed points */ 80 if (numberof(xpp)<numberof(xp)) yp= interp(yp,xpp, xp); 81 82 return yp; 83 } 84