From: dnnant01@ucthpx.uct.ac.za (A.K. Donno) Newsgroups: comp.sources.hp48 Subject: v02i003: leastsq - Least squares curve fitting, Part01/01 Summary: A program that will fit a polynomial to a set a data points Keywords: Least squares Date: 23 Sep 91 03:27:56 GMT Followup-To: comp.sys.hp48 Organization: University of Cape Town Checksum: 2430449214 (verify with brik -cv) Submitted-by: A.K. Donno Posting-number: Volume 2, Issue 3 Archive-name: leastsq/part01 BEGIN_DOC leastsq.doc The feature that I use most in the HP48 stats menu is that of curve fitting, however if you have data that does not fit the standard curve types offered by the calculator, a good fit cannot be obtained. The following program fits a polynomial equation to a set of data points in a similar fashion to the stats menu. It takes a matrix off the stack that has at least two columns in it (this is not checked), allows you to specify which column is x data and which is y data, and then fits a curve to the data using a polynomial of an order that you specify. Please note that more than 2 columns are allowed, but the program only works with the two that are selected. The program uses a standard least squares fit, but operates fairly quickly (4.5 mins to fit a 15th order to 21 data pairs). Note that an option allows you to specify whether the program must plot the fitted curve, or just return the equation. The program does not recalculate the polynomial, unless you change the x or y columns or the order of the polynomial required. The following menu is used: XCOL - takes an integer off the stack and makes that column the x data column YCOL - takes an integer off the stack and makes that column the y data column ORDER - takes an integer off the stack and makes it the order of the poly- nomial that is required PLOT - calculates (if necessary) and plots the fitted polynomial FIT - calculates (if necessary) and returns the polynomial that fits the data EXIT - returns back to the LEAST directory CST menu (you can put whatever you want in this menu although a suggested menu is given) The following is a list of the programs that make up LEAST: LEAST - main program (menus ...) PLOT - plotting function FIT - curve fitter (brains) POLY - converts a list of numbers on the stack into a algebraic polynomial of the order given in level one (the numbers are the coefficients) SETRNG - sets the range for the plot XYR - redraws the top line of the LEAST screen (XCOL: 1 YCOL: 2 ORDER: 3) DELEQ - deletes the current equation EQ The following variables are produced by LEAST: ORDR - order of the polynomial YC - y data column XC - x data column dDAT - summation data (same as the stats menu) EQ - current equation Checksum : # 63743d Bytes : 1950.5 I hope that this program is of use to others, and would appreciatiate it if you would mail any comments to me. END_DOC Here is the ASC file for the program followed by the source: (decode with ASC->) BEGIN_ASC leastsq.asc %%HP: T(3)A(D)F(.); "69A20FF7F390000000303435453047A2084E2050C45414354584E204005C4F44 584E203064944547A20B213047A20C2A20F00000555257454D9D20E163247A20 84E2020541584E20400505142584E2040F425442584E2020953484E202085348 4E204058441445B2130EFE0293632B2130B213047A20C2A20D000054859445D9 D20E1632041A193632B2130B2130B2130F0100504454C4541550D9D20E16323C E224563284E202054159763282EC1683A2D9AE1AFE22D9D204563284E2020541 597632EFE02B21305DF2293632B2130970003085952530D9D20E1632916C11C4 32D6E205066C6167637E16323392010000000000005495D2C133920100000000 00006495D2C13392010000000000007495D2C13392010000000000008495D2C1 3392010000000000009495D2C13392010000000000000595D2C1C2A20F000085 36F6C6A384E20208534B0BC176BA1C2A2011000029536F6C6A376BA184E20209 534B0BC176BA1C2A203100002F427465627A376BA184E2040F4254425B0BC176 BA1C2A2033000A0F5F5F5F5F5F5F5F5F5F5F5F5F5F5F5F5F5F5F5F5F5F576BA1 9C2A2485A19C2A24A5A1D6E205066C6167637F76C1EF53293632B21305D10060 35544525E47460D9D20E1632E7EF184E202085346C7D1E7EF184E202095346C7 D11C432D6E204087D61687D6E204097D61687E16324BEF184E202085346C7D14 BEF184E202095346C7D11C432D6E204087D696E6D6E204097D696E6E1632D6E2 04087D61687D6E204087D696E6126E1D6E204097D61687D6E204097D696E6146 E1EF532EF53293632B2130811004005F4C49540D9D20E16321C432D6E201046E 16320B3A25D2C14563284E20108597632D6E201046D20B1EEDA1D6E2010469C2 A290DA14B2A20A132D6E201046DBBF14563284E20108597632D6E201046D20B1 EEDA176BA1683A208332EF53293632B21308C0003064944530D9D20E163284E2 040584414458B9C1B7FC18DBF18DBF11C432D6E201046E16323CE22D6E201046 84E2040F4254425D5CE1AFE22D9D20D6E20104684E2040F42544259C2A276BA1 ED2A2387C19C2A2681D1ED2A284E2040F42544259C2A276BA10A132D6E201096 9C2A2D6E2010460A132D6E2010A6D6E2010A6D6E201096ED2A2387C184E20405 8441445D6E2010A684E20208534ED2A2387C16C7D1D6E2010969C2A290DA1D20 B1704D1C4232C42329C2A2D6E2010460A132D6E20109684E204058441445D6E2 0109684E20209534ED2A2387C16C7D1C4232D6E2010469C2A2ED2A2387C1900D 1DBBF178BF178BF1293D1DBBF1EEDA1872B1DBBF1293D1EEDA1DBBF1EEDA1B7F C18DBF184E2040F425442584E204005F4C495B21305BF22C2A20D5000542525F 425A302F42746562702E30205F696E64737A0F5F5F5F5F5F5F5F5F5F5F5F5F5F 5F5F5F5F5F5F5F5F55DF224563284E2020541597632DCC02EF53293632B21307 D2004005C4F44540D9D20E1632F52E184E2040584414458B9C1B7FC18DBF18DB F11C432D6E201046E16329C2A2D6E2010460A132D6E20109684E204058441445 D6E20109684E20208534ED2A2387C16C7D184E204058441445D6E20109684E20 209534ED2A2387C16C7D1E97C1A92E1A13E1C42326C1E1091E1EF53293632B21 30EF00050C45414354550D9D20E163284E20504454C45415B2DF116DF1858A19 C2A24563284E2020853497632DCC02ED2A24563284E2020953497632DCC023F2 A24563284E2040F425442597632DCC0284E203085952547A2047A20C2A20D000 08534F4C447A20D9D20E16324563284E2020853497632DCC0284E20504454C45 41584E203085952593632B2130D9D20E16324563284E2020853497632DCC0284 E20504454C4541584E203085952593632B2130D9D20E163284E2020853484E20 3085952593632B2130B2130B213047A20C2A20D00009534F4C447A20D9D20E16 324563284E2020953497632DCC0284E20504454C4541584E203085952593632B 2130D9D20E16324563284E2020953497632DCC0284E20504454C4541584E2030 85952593632B2130D9D20E163284E2020953484E203085952593632B2130B213 0B213047A20C2A20F0000F42544542547A20D9D20E16324563284E2040F42544 2597632DCC0284E20504454C4541584E203085952593632B2130D9D20E163245 63284E2040F425442597632DCC0284E20504454C4541584E203085952593632B 2130D9D20E163284E2040F425442584E203085952593632B2130B2130B213047 A20C2A20D000005C4F445D9D20E16323CE224563284E202054159763282EC168 3A2279E1AFE22D9D2084E206035544525E47484E2030649445B21305DF223CE2 24563284E202054159763282EC1ED2A2279E1AFE22D9D2084E202054159C2A24 85A1ED2A2F17A1B21305BF22D9D20166E184E204005C4F445AB2E1B21305DF22 84E203085952593632B2130B213047A20C2A20B0000649445D9D20E16323CE22 4563284E202054159763282EC1683A2279E1AFE2284E20306494455DF223CE22 4563284E202054159763282EC1ED2A2279E1AFE22D9D2084E202054159C2A248 5A1ED2A2F17A1B21305BF22D9D204563284E202054159763204B02B21305DF22 84E203085952593632B2130B213047A20C2A20D000054859445D9D20E16324B2 A26911293632B2130B2130B2130D511293632B2130FF8F" END_ASC BYTES: #F8FFh 1952.5 BEGIN_UU leastsq.uue begin 644 leastsq.bin M2%!(4#0X+466*O!_/PD````#0U-4`W0J@.0"!4Q%05-42"Y``,7T1(7D`@-&\ M251T*K`2`W0JP*("#P``525U5-39`AXV0J<"2"X@4!2%Y`($4%!!4D@N0/`D" M122%Y`("64-(+B"`-83D`@2%1$%4*S'@[R`Y-K(2`RLQ0*<"+"K0``!%6$E4P MG2W@82-`H9%C(RLQL!(#*S'P$``%1$5,15$%G2W@82/#+D)E(T@N(%`4E6KP:("$0``DC7VQJ9SMAI(+B"0` M-;2P'&>KP:("$P``\B1'5B:G<[8:2"Y`\"1%)+6P'&>KP:(",P"@\/7U]?7U! M]?7U]?7U]?7U]?7U]?7U]76V&LFB0E@:R:)"6AIM+E!@QA9V-O=G'/XUDF,CC M*S%0'0`&4T544DY'!ITMX&$C?OZ!Y`("6$/&U^'G'T@N()`U9'P=P332Y@($M M>&UA>&TN0)#7%H;G82.T_H'D`@)80\;70>L?2"X@D#5D?!W!--+F`@1X;6ENT M;2Y`D->6YN9A(VTN0(#7%H;7Y@($>&UI;B'FT>8"!'EM87AM+D"0UY;F%F0>Z M_C7B7R,Y-K(2`Q@!0`#UQ)1%T-D"'C823"-M+A!`YF$CL*-2+1Q4-H+D`@%8Y M>3;2Y@(!9"VPX=X:;2X00)8L*@FM02LJH#'2Y@(!9+W[064C2"X0@)5G(VTNI M$$#6`AONK7&V&H:C`C@C_C628R,K,8`,``-&250#G2W@82-(+D!02!1$A9L\^!O1_8^Q%,(VTN$$#F82/#+M+F`@%D2"Y`\"1%)-7%'OHNTMD";2X00(;DA M`@1/4D12R:)RMAK>HC)X',FB8A@=WJ*"Y`($3U)$4LFB!Q(+D!02!1$U>8"`6I(+B"`->0M, M*H/'87P=;2X0D)8L*@FMT0(;!]3!)"-,,I(L*FTN$$`&&B-M+A"0AN0"!(5$P M051M+A"0AN0"`EE#WJ(R>!S&U\$D(VTN$$"6+"K>HC)X'`G0T;L?A_MQN!^22 MT]&['^ZM@2<;O?LA.1WNK=&['^ZML?<!Y`($A41!5&TN$)"&Y`("64/>HC)X',;7M MX7D-H+D`@5$14Q%W M42O]$=8?6*B1+"I4-H+D`@)80WDVTLP@WJ)"92-(+B"0-91G(\T,,B\J5#:"2 MY`($3U)$4GDVTLP@2"XP@)4E1:<"="K`H@(-`(`U],1$IP*=+>!A(U0V@N0"V M`EA#>3;2S"!(+E!`5,14%(7D`@-865(Y-K(2`YTMX&$C5#:"Y`("6$-Y-M+,. M($@N4$!4Q%04A>0"`UA94CDVLA(#G2W@82-(+B"`-83D`@-865(Y-K(2`RLQV ML!(#="K`H@(-`)`U],1$IP*=+>!A(U0V@N0"`EE#>3;2S"!(+E!`5,14%(7D1 M`@-865(Y-K(2`YTMX&$C5#:"Y`("64-Y-M+,($@N4$!4Q%04A>0"`UA94CDVL MLA(#G2W@82-(+B"0-83D`@-865(Y-K(2`RLQL!(#="K`H@(/`/`D150D1:<"+ MG2W@82-4-H+D`@1/4D12>3;2S"!(+E!`5,14%(7D`@-865(Y-K(2`YTMX&$CW M5#:"Y`($3U)$4GDVTLP@2"Y00%3$5!2%Y`(#6%E2.3:R$@.=+>!A(T@N0/`DS M122%Y`(#6%E2.3:R$@,K,;`2`W0JP*("#0``Q?1$U=D"'C8R[")4-H+D`@)%6 M47DV@N(\BG2V`Y`("15')HD)8&MZB\G$:*S%0^R*=+1!F'D@N@ M0`#%]$2E*QXK,5#](D@N,("5)95C(RLQL!(#="K`H@(+`&"41-79`AXV,NPB' M5#:"Y`("15%Y-H+B'(:C(I<>^BZ"Y`(#1DE4U2\R[")4-H+D`@)%47DV@N(<: MWJ(BEQ[Z+M+9`D@N(%`4E2PJA*7A+2H?I[$2`[4OTMD"5#:"Y`("15%Y-@*TB M("LQ4/TB2"XP@)4EE6,C*S&P$@-T*L"B`@T`4(251-79`AXV0BLJEA&28R,KU .,;`2`RLQT!4A.3:R$@.R# `` end END_UU BEGIN_RPL leastsq.rpl %%HP: T(3)A(D)F(.); DIR LEAST \<< DELEQ CL\GS \GS+ CLLCD 1 'XC' STO 2 'YC' STO 3 'ORDR' STO XYR { { "XCOL" { \<< 'XC' STO DELEQ XYR \>> \<< 'XC' STO DELEQ XYR \>> \<< XC XYR \>> } } { "YCOL" { \<< 'YC' STO DELEQ XYR \>> \<< 'YC' STO DELEQ XYR \>> \<< YC XYR \>> } } { "ORDER" { \<< 'ORDR' STO DELEQ XYR \>> \<< 'ORDR' STO DELEQ XYR \>> \<< ORDR XYR \>> } } { "PLOT" \<< IF 'EQ' VTYPE -1 == THEN SETRNG FIT END IF 'EQ' VTYPE 2 == THEN EQ 1 DISP 2 WAIT ELSE FUNCTION PLOT GRAPH END XYR \>> } { "FIT" \<< IF 'EQ' VTYPE -1 == THEN FIT END IF 'EQ' VTYPE 2 == THEN EQ 1 DISP 2 WAIT ELSE 'EQ' RCL END XYR \>> } { "EXIT" \<< 0 MENU \>> } } TMENU \>> PLOT \<< ERASE \GSDAT SIZE OBJ\-> DROP DROP \-> d \<< 1 d FOR i \GSDAT i XC 2 \->LIST GET \GSDAT i YC 2 \->LIST GET R\->C C\->PX PIXON NEXT DRAX DRAW \>> \>> FIT \<< \GSDAT SIZE OBJ\-> DROP DROP \-> d \<< IF d ORDR > THEN d ORDR 1 + 2 \->LIST 1 CON 2 ORDR 1 + FOR i 1 d FOR j j i 2 \->LIST \GSDAT j XC 2 \->LIST GET i 1 - ^ PUT NEXT NEXT 1 d FOR i \GSDAT i YC 2 \->LIST GET NEXT d 1 2 \->LIST \->ARRY SWAP DUP DUP TRN SWAP * INV SWAP TRN * SWAP * OBJ\-> DROP ORDR POLY ELSE "ERROR: Order > Points ______________________" END 'EQ' STO \>> \>> POLY \<< \-> d \<< -3 CF 'X' d ^ * d 1 - 0 FOR d SWAP 'X' d ^ * + -1 STEP \>> \>> SETRNG \<< MAX\GS XC GET MAX\GS YC GET \-> xmax ymax \<< MIN\GS XC GET MIN\GS YC GET \-> xmin ymin \<< xmax xmin XRNG ymax ymin YRNG \>> \>> \>> XYR \<< RCLF \-> flags \<< -45 CF -46 CF -47 CF -48 CF -49 CF -50 CF "Xcol:" XC \->STR + " Ycol:" + YC \->STR + " Order:" + ORDR \->STR + " ______________________" + 1 DISP 1 FREEZE flags STOF \>> \>> DELEQ \<< IF 'EQ' VTYPE -1 \=/ THEN 'EQ' PURGE END \>> CST { LEAST PLOT FIT { } { "PURGE" \<< { EQ PPAR ORDR YC XC \GSDAT } PURGE \>> } { "EXIT" \<< HOME \>> } } END END_RPL ============================================================================== Anthony Donno | Department of Electrical Engineering | Caesar entered on his head University of Cape Town | he had a helmet on each foot New South Africa | he had a sandal in his hand | he had his trusty sword to boot e-mail : dnnant01@ucthpx.uct.ac.za | phone : +27 +21 6892420 | ==============================================================================