DIRECTORYLeastSq, Random, Real;LeastSqImpl: CEDAR PROGRAM IMPORTS Random, Real EXPORTS LeastSq =BEGINNewPointSet: PUBLIC PROC [ dim: NAT ] RETURNS [ ps: LeastSq.PointSet ] = {ps _ NEW[LeastSq.PointSetRec[dim]];ps.count _ 0;FOR row: NAT IN [0..dim) DOps[row] _ NEW[LeastSq.RowRec[dim]];FOR col: NAT IN [0..dim) DOps[row][col] _ 0.0;ENDLOOP;ENDLOOP;};AddPoint: PUBLIC PROC [ ps: LeastSq.PointSet, pt: LeastSq.Row ] = {IF pt.colDim # ps.rowDim THEN ERROR;ps.count _ ps.count+1;FOR row: NAT IN [0..ps.rowDim) DOFOR col: NAT IN [0..row] DOps[row][col] _ ps[row][col]+pt[row]*pt[col];ENDLOOP;ENDLOOP;};LeastSqFit: PUBLIC PROC [ ps: LeastSq.PointSet ] RETURNS [ eq: LeastSq.Row ] = {basis: LeastSq.PointSet = NEW[LeastSq.PointSetRec[ps.rowDim]];FOR row: NAT IN [0..ps.rowDim) DO -- fill ps symmetricallyFOR col: NAT IN (row..ps.rowDim) DOps[row][col] _ ps[col][row];ENDLOOP;ENDLOOP;FOR row: NAT IN [0..ps.rowDim) DOinnerProd: REAL;col: NAT;basis[row] _ NEW[LeastSq.RowRec[ps.rowDim]];FOR col IN [0..ps.rowDim) DO basis[row][col] _ 0.0 ENDLOOP;FOR prevrow: NAT IN [0..row) DOinnerProd _ 0.0;FOR col IN [0..ps.rowDim) DOinnerProd _ innerProd+basis[prevrow][col]*ps[row][col];ENDLOOP; -- colFOR col IN [0..ps.rowDim) DObasis[row][col] _ basis[row][col]-innerProd*basis[prevrow][col];ENDLOOP; -- colENDLOOP; -- prevrowbasis[row][row] _ 1.0;IF row<ps.rowDim-1 THEN -- normalize this new orthogonal polynomialBEGINinnerProd _ 0.0;FOR col IN [0..ps.rowDim) DOFOR col1: NAT IN [0..ps.rowDim) DOinnerProd _ innerProd+basis[row][col]*basis[row][col1]*ps[col][col1];ENDLOOP; -- col1ENDLOOP; -- colIF innerProd>0.0 THEN {normFact: REAL = 1.0/Real.SqRt[innerProd];FOR col IN [0..ps.rowDim) DObasis[row][col] _ normFact*basis[row][col];ENDLOOP; -- col};END;ENDLOOP; -- roweq _ basis[ps.rowDim-1];};NewRow: PROC [ dim: NAT ] RETURNS [ pt: LeastSq.Row ] = {pt _ NEW[LeastSq.RowRec[dim]]};Test: PROC [ eqIn: LeastSq.Row, n: INT, noise: REAL _ 0 ] RETURNS [ eqOut: LeastSq.Row ] =BEGINps: LeastSq.PointSet = NewPointSet[eqIn.colDim];pt: LeastSq.Row = NEW[LeastSq.RowRec[eqIn.colDim]];pt[0] _ 1.0;FOR i: INT IN [0..n) DOsum: REAL _ eqIn[0]+noise*Random.Choose[-10000, 10000];FOR j: NAT IN [1..eqIn.colDim-1) DOpt[j] _ Random.Choose[-10000, 10000];sum _ sum+eqIn[j]*pt[j];ENDLOOP;pt[eqIn.colDim-1] _ -sum;AddPoint[ps, pt];ENDLOOP;eqOut _ LeastSqFit[ps];END;END.  °  LeastSqImpl.mesaLast Edited by: McCreight, April 12, 1985 11:49:30 am PSTAt this point, ps[row][col] contains the sum over all input points of pt[row]*pt[col], where pt[k] =1 if k=0,an independent variable if k>0 and k<dim-1,the dependent variable if k=dim-1.We begin with a basis set of polynomials: 1, ind1, ind2, ....  We use the Gram-Schmidt algorithm (see the textbook Numerical Analysis, second edition, by Johnson and Riess, Addison Wesley, 1982, pages 255 and following for details) to construct a new basis set of polynomials that is orthogonal with respect to the given set of independent variable pairs, and normal in magnitude (in that the Euclidean norm of each polynomial is one).The best least square approximation (with respect to the original basis set) is then a linear combination of the orthonormal polynomials, where the coefficient of polynomial i is the inner product of the dependent values with the vector of values of that polynomial evaluated at the independent variables. Э╩H  Ш JЪЬЩJЩ9JШ Ъ╧k	Ш	JЪЬШЧJШ Ъ	ЬЭЬЭЬЭЬЭЬ
ШAJЪЭШJШ Ъ
╧nЬЭЬЭЬЭЬЭЬШJJЪЬЭЬШ#JЪЬШЪЭЬЭЬЭЬ
ЭШJЪЬ
ЭЬШ#ЪЭЬЭЬЭЬ
ЭШJЪЬШJЪЭЬШЧJЪЭЬШЧJЪЬШJШ ЧЪЮЬЭЬЭЬ.ШCJЪЭЬЭЬЭЬШ$JЪЬШЪЭЬЭЬЭЬЭШ!ЪЭЬЭЬЭЬ
ЭШJЪЬ,Ш,JЪЭЬШЧJЪЭЬШЧJЪЬШJШ ЧЪЮ
ЬЭЬЭЬШPЩdJЩ	JЪЬ+Щ+JЪЬ"Щ"JЩ ЧJЩ│JЩ JЩ▒JШ JЪЬЭЬ!Ш>JШ Ъ	ЭЬЭЬЭЬЭЬ╧cШ:ЪЭЬЭЬЭЬЭШ#JШJЪЭЬШЧJЪЭЬШJШ ЧЪЭЬЭЬЭЬЭШ!JЪЬЭЬШJЪЬЭЬШ	JЪЬЭЬШ,JЪЭЬЭЬЭЬЭЬШ;ЪЭЬ
ЭЬЭЬ
ЭШJЪЬШЪЭЬЭЬЭШJЪЬ7Ш7JЪЭЬЯШЧЪЭЬЭЬЭШJЪЬ@Ш@JЪЭЬЯШЧJЪЭЬЯ
ШJШ ЧJШJШ ЪЭЬЭЬЯ+ШCJЪЭШJЪЬШЪЭЬЭЬЭШЪЭЬЭЬЭЬЭШ"JЪЬEШEJЪЭЬЯШЧJЪЭЬЯШЧЪЭЬЭЬШJЪЬ
ЭЬШ*ЪЭЬЭЬЭШJЪЬ+Ш+JЪЭЬЯШЧJШЧJЪЭЬШЧJЪЭЬЯШJШ JЪЬШЧJШJШ JШ ЧЪЮЬЭЬЭЬЭЬШ9JЪЬЭЬШJЪЬШJШ ЧЪ
ЮЬЭЬЭЬ	ЭЬЭЬШZJЪЭШJЪЬ0Ш0JЪЬЭЬШ3JШЪЭЬЭЬЭЬЭШJЪЬЭЬ.Ш7ЪЭЬЭЬЭЬЭШ#JЪЬ%Ш%JЪЬШJЪЭЬШЧJЪЬШJЪЬШJЪЭЬШЧJЪЬШJЪЭЬШJШ ЧЪЭЬШШ JШ JШ ЧЧЧЧ ЕЧ    	Д  ─  