# hp41programs

least-squares Least-Squares Approximation for the HP-41

Overview

1°)  The Least-squares Line(s)

a)  Linear Regression
b)  Another Least-squares Line
c)  Both Least-squares Lines

2°)  Least-Squares Parabola
3°)  Least-Squares Polynomials: equally-spaced abscissas
4°)  Linear Combination of 3 Functions
5°)  Linear Combination of n Functions
6°)  Least-Squares Plane
7°)  Linear Combination of n Functions ( Generalization to a Surface )
8°)  Least-Squares 3D-Hyperplane

-Given a set of N data points ( xi , yi ) , we seek a best fitting curve y = a1 f1(x) + a2 f2(x) + ...... + an fn(x)           ( n <= N )
where  f1 ,  f2 , ...... ,  fn      are n  known functions and  a1 , a2 , ...... , an   are unknown.
-The least-squares principle minimizes the sum of "vertical distances" to a curve:        SUMi   [ yi - a1 f1(xi) + a2 f2(xi) + ...... + an fn(xi) ]2
leading to a n*n linear system called the "normal equations".

1°)  The Least-Squares Line(s)

a)  Linear Regression  ( SigmaREG 00 )

-Here, the curve is a straight line y = m.x + p  given by:

m = ( Sum xi yi - n µx µy ) / ( Sum xi2 - ( Sum xi )2 / n )                r = the linear correlation coefficient = ( Sum xi yi - n µx µy ) / (( n-1) sx sy)
p =   µy - m.µx                                                                and       µx , µy ,  sx , sy    are the means and standard deviations of x and y.

Data Registers:                                               ( Registers R00 thru R05 are to be initialized before executing "LR" )

•   R00 = Summation of x values                •   R02 = Summation of y values                  •   R04 = Summation of x.y values
•   R01 = Summation of x2 values               •   R03 = Summation of y2 values                •   R05 = Number of data points

Flags: /
Subroutines: /

 01  LBL "LR" 02  STO M   03  MEAN 04  * 05  RCL 05   06  * 07  RCL 04 08  - 09  STO Z 10  SDEV 11  RCL 05   12  DSE X 13  * 14  * 15  / 16  CHS 17  X<>Y 18  RCL 00   19  X^2 20  RCL 05 21  / 22  RCL 01 23  - 24  / 25  MEAN 26  LASTX 27  * 28  ST- Y 29  X<> L 30  RCL M   31  SIGN 32  RDN 33  ST* M 34  X<>Y 35  ST+ M 36  X<>Y 37  0 38  X<> M   39  END

( 55 bytes / SIZE 006 )

 STACK INPUTS OUTPUTS T / r Z / p Y / m X x y(x) = m.x + p L / x

Example:     Find the least-squares linear function for the following data:

 y 2 3 3 5 5 x 1 2 3 4 5

and evaluate the linear estimates   y(2.5)  and  y(7)

a-  [SREG 00]   [CLS]        ( SigmaREG 00 and CLSigma )

b-  2  ENTER^     1      S+            (  Sigma+ )               ( enter y-values first )
3  ENTER^      2     S+
3  ENTER^             S+
5  ENTER^      4     S+
5  ENTER^             S+

c- 2.5 XEQ "LR"  yields   3.2                thus,        y(2.5) = 3.2
RDN              0.8
RDN              1.2               whence     y = 0.8 x + 1.2
RDN              0.9428         whence      r = 0.9428

7     R/S   >>>    y(7) = 6.8

Notes:

-If your favorite statistics registers are, for instance, R11 thru R16 , replace R00 thru R05 by these registers in the above listing.
-Synthetic register M can be replace by any other unused data register.
-This program calculates m , p , r each time: thus, no other data register is used.

b)  Another Least-Squares Line  ( SigmaREG 00 )

-The following program minimizes the sum of squares of perpendicular distances instead of "vertical distances" to a line y = m.x + p
this leads to a quadratic equation and m and p are determined by:

m = -A/2 + sign (r) ( 1 + A2/4 )1/2             where  A = ( sx/sy - sy/sx ) / r    ;    r =  linear correlation coefficient  = ( Sum xi yi - n µx µy ) / (( n-1) sx sy)
p =   µy - m.µx                                                                                    and       µx , µy ,  sx , sy    are the means and standard deviations of x and y.

Data Registers:                                                ( Registers R00 thru R05 are to be initialized before executing "LSL" )

•   R00 = Summation of x values                 •   R02 = Summation of y values                 •   R04 = Summation of x.y values
•   R01 = Summation of x2 values               •   R03 = Summation of y2 values                •   R05 = Number of data points

Flags: /
Subroutines: /

 01  LBL "LSL" 02  STO M 03  SDEV 04  * 05  RCL 05     06  DSE X 07  * 08  RCL 04 09  RCL 00 10  RCL 02 11  * 12  RCL 05     13  / 14  - 15  X<>Y 16  / 17  RDN 18  SDEV 19  / 20  ENTER^ 21  1/X 22  - 23  R^ 24  ST+ X 25  / 26  STO Y       27  X^2 28  1 29  + 30  SQRT 31  R^ 32  SIGN 33  * 34  + 35  MEAN 36  LASTX      37  * 38  ST- Y 39  X<> L 40  RCL M 41  SIGN 42  RDN 43  ST* M 44  X<>Y 45  ST+ M 46  X<>Y 47  0 48  X<> M      49  END

( 67 bytes / SIZE 006 )

 STACK INPUTS OUTPUTS T / r Z / p Y / m X x y(x) = m.x + p L / x

Example:    Find this second least-squares line with the same data:

 y 2 3 3 5 5 x 1 2 3 4 5

and evaluate the linear estimates   y(2.5)  and  y(7)    ( skip to step c if registers R00 thru R05 are unchanged )

a-  [SREG 00]   [CLS]        ( SigmaREG 00 and CLSigma )

b-  2  ENTER^     1      S+            ( Sigma+ )               ( enter y-values first )
3  ENTER^      2     S+
3  ENTER^             S+
5  ENTER^      4     S+
5  ENTER^             S+

c- 2.5 XEQ "LSL"  yields   3.1799                thus          y(2.5) = 3.1799
RDN              0.8402
RDN              1.0794               whence     y = 0.8402 x + 1.0794
RDN              0.9428               whence      r = 0.9428

7     R/S   -->     y(7) = 6.9608

Notes:

-If your favorite statistics registers are, for instance, R11 thru R16 , replace R00 thru R05 by these registers in the above listing.
-Synthetic register M can be replace by any other unused data register.
-This program calculates m , p , r each time: thus, no other data register is used.

c)  Both Lines  ( SigmaREG 00 )

-The 2 programs above may be combined in a single routine.

Data Registers:

R00 = Summation of x values            R02 = Summation of y values         R04 = Summation of x.y values     R06 = m
R01 = Summation of x2 values           R03 = Summation of y2 values       R05 = Number of data points        R07 = p

Flag:   CF 01 for the "classical" linear regression
SF 01 for the second least-squares line ( LSL )

Subroutines: /

 01  LBL A 02  ST+ 05 03  STO 06 04  X<>Y 05  * 06  STO 07   07  ST+ 02 08  LASTX 09  * 10  ST+ 03 11  X<>Y 12  ENTER^ 13  ST* 07 14  RCL 06 15  * 16  ST+ 00 17  * 18  ST+ 01 19  RCL 07 20  ST+ 04 21  RCL 06 22  RTN 23  LBL "LR" 24  RCL 04   25  RCL 00 26  RCL 02 27  * 28  RCL 05 29  / 30  - 31  RCL 05 32  1 33  - 34  / 35  STO 06   36  STO T 37  RDN 38  SDEV 39  ST/ Z 40  X^2 41  ST/ 06 42  RDN 43  ST/ Y 44  FC? 01 45  GTO 01 46  X^2 47  R^ 48  - 49  R^ 50  STO 06 51  ST+ X 52  / 53  STO 07   54  X^2 55  1 56  + 57  SQRT 58  RCL 06 59  SIGN 60  * 61  RCL 07 62  + 63  STO 06 64  LBL 01 65  CLX 66  RCL 00 67  RCL 06 68  * 69  RCL 02 70  X<>Y 71  - 72  RCL 05   73  / 74  STO 07 75  RCL 06 76  RTN 77  LBL F 78  RCL 06 79  * 80  RCL 07 81  + 82  RTN 83  GTO F 84  END

( 109 bytes / SIZE 008 )

 STACK INPUTS OUTPUTS T / cov(x,y) Z / r Y / p X / m

where the least-squares line is   y = mx + p  ,   r = coefficient of correlation  ,  cov(x,y) = sample covariance

Example:

 xi 2 5 10 16 21 yi 4 7 12 19 24 ni 3 4 7 5 2

-Clear registers  R00 thru R05  ( for instance execute CLS )
-Then,  xi  ENTER^   yi  ENTER^    ni   S+  ( [A]  in user mode )  for all the points

2  ENTER^      5  ENTER^      10  ENTER^     16  ENTER^    21  ENTER^
4  ENTER^      7  ENTER^      12  ENTER^     19  ENTER^    24  ENTER^
3  [A]               4  [A]                7   [A]               5   [A]              2   [A]              in user mode

-Then, simply press  R/S  ( or XEQ "LR" )  it yields:

CF 01                                                           SF 01

m  = 1.0694                                                    m  =  1.0793
RDN            p  =  1.6130                               RDN             p  =  1.6039
RDN            r   =  0.9992                               RDN             r   =  0.9992
RDN   cov(x,y) = 38.0143                              RDN   cov(x,y) = 38.0143

-Finally, for a linear estimate - for instance y(7)

7   R/S  ( or [F]  in user mode )  gives  y(7) = 9.0987    if  CF 01  or  y(7) = 9.0958    if  SF 01

Note:

-If all the ni's = 1 , lines 01 thru 22 are unuseful and the data may be stored by   yi  ENTER^  xi   S+   for all the points

2°)  Least-Squares Parabola

-Given a set of n data points ( xi , yi ) , we seek the best fitting parabola  y = a x2 + b x + c

-We have to solve the 3x3 linear system

a Sum x4 + b Sum x3 + c Sum x2 = Sum x2y
a Sum x3 + b Sum x2 + c Sum x  = Sum x.y
a Sum x2 + b Sum x  + c n          = Sum y

-And the correlation coefficient is given by   r = [ a Sum x2y + b Sum x.y + c Sum y - (Sum y)2/n ] / [ Sum y2 - (Sum y)2/n ]

Data Registers:   R00 = Sum x     R02 = Sum y      R04 = Sum x.y    R06 = Sum x3     R08 = Sum x2y        R12 = r
R01 = Sum x2    R03 = Sum y2    R05 = n               R07 = Sum x4     R09 = a , R10 = b , R11 = c
Flags: /
Subroutines: /

 01  LBL I   02  SREG 03    03  CLS         04  SREG 00   05  CLS       06  RTN   07  LBL A   08  S+        09  LASTX   10  X^2   11  ST* L   12  ST* Y   13  ST* X   14  ST+ 07   15  RDN   16  ST+ 08   17  LASTX   18  ST+ 06   19  RCL 05   20  RTN   21  LBL "LSP"   22  RCL 05   23  RCL 07   24  *   25  RCL 01 26  X^2   27  -   28  STO 09   29  RCL 04   30  RCL 05        31  *   32  RCL 00   33  RCL 02   34  *   35  -   36  STO 11   37  *   38  STO 10   39  RCL 05   40  RCL 06   41  *   42  RCL 00   43  RCL 01   44  *   45  -   46  STO 12   47  RCL 05   48  RCL 08   49  *   50  RCL 01 51  RCL 02   52  *   53  -   54  *   55  ST- 10   56  LASTX   57  RCL 01   58  RCL 05         59  *   60  RCL 00   61  X^2   62  -   63  ST* 09   64  *   65  RCL 11   66  RCL 12   67  *   68  -   69  RCL 09   70  RCL 12   71  X^2   72  -   73  ST/ 10   74  /   75  STO 09 76  RCL 01   77  *   78  RCL 00   79  RCL 10   80  *   81  +   82  RCL 02   83  X<>Y   84  -   85  RCL 05         86  /   87  STO 11   88  RCL 02   89  *   90  RCL 04   91  RCL 10   92  *   93  +   94  RCL 08   95  RCL 09   96  *   97  +   98  RCL 03   99  RCL 02 100  X^2 101  RCL 05 102  / 103  ST- Z 104  - 105  / 106  STO 12 107  RCL 11 108  RCL 10       109  RCL 09 110  RTN 111  LBL F 112  STO Y 113  RCL 09 114  * 115  RCL 10 116  + 117  * 118  RCL 11 119  + 120  RTN 121  GTO F 122  END

( 148 bytes / SIZE 013 )

 STACK INPUTS#i OUTPUTS#i OUTPUTS T / / r Z / / c Y yi / b X xi i a

where   y = a x2 + b x + c  is the best fitting parabola and  r = coefficient of correlation

Example:    Given the following data points

 yi 0 3 4 5 1 xi 0 1 2 3 5

-First, we clear registers R00 thru R08   XEQ I  or  press  [I]  in user mode or  RTN  R/S
-We store the data

0   ENTER^        3   ENTER^       4   ENTER^       5   ENTER^           1   ENTER^
0   [A]                 1   [A]                2   [A]                3   [A]                    5   [A]              all in user mode

-Then press R/S  ( or XEQ "LSP" )  it yields   a = -0.6701  = R09
RDN     b =  3.5670  = R10
RDN     c = -0.0206  = R11
RDN     r =   0.9808  = R12

So the least-squares parabola is  y = -0.6701 x2 + 3.5670 x - 0.0206

-To obtain an estimation of, say  y(4) & y(7)   press

4  R/S  ( or [F]  in user mode )   y(4) =  3.5258
7  R/S  ( or [F]  in user mode )   y(7) = -7.8866   ... and so on ...

3°)  Least-Squares Polynomials , Equally-Spaced Abscissas

-Given a set of (N+1) data points ( x0 , y0 ) ,  ( x1 , y1 )  , ........................ ,  ( xN , yN ) ,
you can use the program of paragraph 5 to find the least-squares polynomial p(x) of degree m  ( m <= N )

-However, the normal equations generally involve a very ill-conditioned matrix for large m- and N-values.
-So, another approach is used, but we assume that the xi's are equally spaced
-We define orthogonal polynomials Pk,N(t)  by                                                                       ( t = (x-x0)/h  with  h = xi+1 - xi )

Pk,N(t) = SUM i=0,1,...,k (-1)i  Cki Ck+ii  t(i)/N(i)

where   t(i) = t(t-1)(t-2)....(t-i+1)  and  N(i) = N(N-1)(N-2).....(N-i+1)      ( = 1 if  i = 0 )
and    Cki = k!/(i!(k-i)!)

-The polynomial  p(x) can then be written  p(x) = a0 P0,N(t) + a1 P1,N(t) + .................... + am Pm,N(t)

with   ak = [ SUM t=0,1,...,N  yt Pk,N(t) ]/[ SUM t=0,1,...,N  (Pk,N(t))2 ]

-These coefficients are computed and stored when the program reaches line 76
-Then ( lines 76 to 193 ), the program calculates the coefficients in the more standard form:   c0 + c1 t + c2 t2 + ...... + cm tm

-In fact, the Pk,N(t) are computed by the recurrence relations:

P-1,N(t) = 0 ,  P0,N(t) = 1 ,  (k+1)(N-k) Pk+1,N(t) = (2k+1)(N-2t) Pk,N(t) - k(N+k+1) Pk-1,N(t)

Data Registers:           •  R00 = N  ( there are N+1 data points )     ( Registers R00 and R10 thru R10+N are to be initialized before executing "LSPN" )

R01 thru R09: temp

•  R10 = y0  •  R11 = y1  •  R12 = y2   .........................   •  R10+N = yN

>>>>    at the end          R11+N = a0  ,  R12+N =  a1  ,  .............................. ,  R11+N+m = am
R12+N+m = c0  , R13+N+m = c1  , ........................... ,  R12+N+2m = cm

R13+N+2m to R13+N+3m = the coefficients of  Pm-1,N (t)
R14+N+3m to R14+N+4m = the coefficients of  Pm,N (t)

Flags: /
Subroutines: /

-If you don't have an HP-41CX, replace line 89 ( CLRGX ) by

STO 02           LBL 00                GTO 00
SIGN             STO IND 02         X<> L
CLX               ISG 02

 01  LBL "LSPN"   02   E3   03  /   04  STO 01   05  RCL 00   06  11   07  +   08  STO 09   09  LBL 01   10  RCL 00   11  STO 02   12  SIGN   13  STO 05   14  RCL 10   15  STO 04   16  LBL 02   17  RCL 01   18  INT   19   E3   20  /   21  STO 03   22  CLST   23  SIGN   24  LBL 03   25  ISG 03   26  FS? 30   27  GTO 04   28  X<>Y   29  RCL 00   30  RCL 03   31  INT   32  STO 06   33  STO 07   34  ST+ Y   35  DSE X   36  ST+ 07   37  STO 08   38  *   39  * 40  RCL 00   41  ST- 08   42  RCL 02   43  ST+ X   44  -   45  RCL 07   46  *   47  R^   48  *   49  -   50  RCL 06   51  RCL 08   52  *   53  /   54  GTO 03   55  LBL 04   56  X^2   57  ST+ 05   58  RCL 02   59  10   60  ST+ Y   61  RCL IND Y    62  LASTX   63  *   64  ST+ 04   65  DSE 02   66  GTO 02   67  RCL 01   68  RCL 09   69  +   70  RCL 04   71  RCL 05   72  /   73  STO IND Y   74  ISG 01   75  GTO 01   76  RCL 01   77  INT   78  ENTER^ 79  STO Z   80  4   81  *   82  DSE X   83  RCL 09   84  ST+ Z   85  +   86   E3   87  /           88  +   89  CLRGX      90  INT         91  RCL IND 09   92  STO IND Y   93  RDN   94  +   95  STO 04   96  +   97  ENTER^   98  ENTER^   99  SIGN 100  STO IND Y 101  + 102   E3 103  / 104  + 105  STO 05 106  1.999 107  STO 03 108  RCL 01 109  FRC 110  ISG X 111  STO 01 112  CLX 113  STO 02 114  LBL 05 115  RCL 01 116  INT 117  STO 06 118  1 119  - 120  STO 07 121  RCL 00 122  RCL 06 123  + 124  * 125  RCL IND 04 126  * 127  RCL 06 128  RCL 07 129  + 130  RCL IND 05 131  * 132  STO Z 133  RCL 00 134  ST- 07 135  * 136  - 137  RCL 06 138  RCL 07 139  * 140  ST/ Z 141  / 142  X<>Y 143  ST+ X 144  X<> 02 145  + 146  X<> IND 05 147  STO IND 04 148  ISG 04 149  CLX 150  ISG 05 151  GTO 05 152  RCL 03 153  ISG X 154  INT 155  ST- 04 156  ST- 05 157  RCL 01 158  RCL 04 159  RCL 09 160  ST+ Z 161  + 162  2 163  / 164  RCL IND Y 165  SIGN 166  LBL 06 167  CLX 168  RCL IND 05 169  LASTX 170  * 171  ST+ IND Y 172  ISG Y 173  CLX 174  ISG 05 175  GTO 06 176  RCL 03 177  ST- 05 178  ISG 03 179  ISG 01 180  GTO 05 181  INT 182  RCL 09 183  + 184   E3 185  / 186  RCL 09 187  + 188  STO Y 189  RCL 01 190   E-3 191  + 192  + 193  END

( 268 bytes / SIZE N+4m+15 )

 STACK INPUTS OUTPUTS Y / bbb.eee(ai) X m bbb.eee(ci)

Example:

 t 0 1 2 3 4 5 6 y 1.2 3.4 5.1 6.7 7.6 8 8.4

N = 6  STO 00

1.2  STO 10         6.7  STO 13        8.4  STO 16
3.4  STO 11         7.6  STO 14
5.1  STO 12           8   STO 15

-If we seek a polynomial of degree 4

4  XEQ "LSPN"  >>>>  22.026 = control number of the coefficients  ci's            --- execution time = 1mn50s ---
RDN   17.021 = control number of the coefficients  ai's

and we get:

p(t) = 5.771429 - 3.567857 P1,N(t) - 1.005952 P2,N(t) - 0.016667 P3,N(t) + 0.037013 P4,N(t)
p(t) = 1.217965 + 2.088023 t + 0.093561 t2 - 0.083586 t3 + 0.007197 t4

Notes:

-Unlike "LSF" ( cf §5 ),  this routine can find least-squares polynomials of degree m = 20 or more - provided N+4m+15 doesn't exceed 319
-With m = N , you have the Lagrange Polynomial.

-Since the Pk,N(t) are orthogonal, we obtain least-squares approximations of lower degree by truncation.
-So, you don't have to execute the whole program again to find the polynomial of degree 3:

p(t) = 5.771429 - 3.567857 P1,N(t) - 1.005952 P2,N(t) - 0.016667 P3,N(t)

-And to get the standard expression, simply add

GTO 12       ISG X                  after line 75
LBL 10          +
.1               STO 01
%              LBL 12

then   3  XEQ 10  gives:   p(t) = 1.180952 + 2.451984 t - 0.226190 t2 + 0.002778 t3   in 19 seconds

-However, the coefficient a4 is lost!

-This program has been used to find 3 polynomials of degree 6 that best fit the coordinates of Pluto from 62 positions between 2005 and 2010
( cf "Astronomical Ephemeris for the HP-41" )
-The execution time was about 30 minutes/coordinate with a "true" HP-41, but less than 3 seconds with a V41 emulator!

4°)  Linear Combination of 3 functions

-Given a set of data points ( xi , yi ), we seek a best fitting curve y = a f(x) + b g(x) + c h(x)
where f , g , h are 3 known functions and a , b , c are unknown.

Data Registers:           ( Registers R01 R02 R03 are to be initialized and R07 thru R15 are to be cleared before executing this program )

•  R01 = f name       R04 = a      R07 = Sum f.f       R10 = Sum g.g      R13 = Sum y.f
•  R02 = g name      R05 = b      R08 = Sum f.g      R11 = Sum g.h      R14 = Sum y.g              ( R00 and R16: temporary data storage )
•  R03 = h name      R06 = c      R09 = Sum f.h      R12 = Sum h.h       R15 = Sum y.h

Flags:  /
Subroutines:   the 3 functions f , g , h  ( assuming x is in X-register upon entry )

 01  LBL A   02  X<>Y   03  STO 00   04  X<>Y   05  STO 16   06  XEQ IND 01   07  STO 04   08  X^2   09  ST+ 07   10  LASTX   11  RCL 00   12  *   13  ST+ 13   14  RCL 16   15  XEQ IND 02   16  STO 05   17  X^2   18  ST+ 10   19  LASTX   20  RCL 00   21  *   22  ST+ 14   23  RCL 16   24  XEQ IND 03   25  STO 06   26  X^2   27  ST+ 12   28  LASTX   29  RCL 00   30  *   31  ST+ 15   32  RCL 04 33  RCL 05   34  *   35  ST+ 08   36  RCL 06   37  LASTX   38  *   39  ST+ 11   40  RCL 04   41  RCL 06   42  *   43  ST+ 09   44  RCL 16   45  RTN   46  LBL "LSF3"    47  RCL 07    48  RCL 10   49  RCL 12   50  *   51  RCL 11   52  X^2   53  -   54  STO 04   55  *   56  RCL 09   57  RCL 11   58  *   59  RCL 08   60  RCL 12   61  *   62  -   63  STO 05   64  RCL 08 65  *   66  +   67  RCL 08   68  RCL 11   69  *   70  RCL 09   71  RCL 10   72  *   73  -   74  STO 06   75  RCL 09   76  *   77  +   78  STO 00           79  RCL 04   80  RCL 13   81  *   82  RCL 05   83  RCL 14   84  *   85  +   86  RCL 06   87  RCL 15   88  *   89  +   90  RCL 00   91  /   92  STO 04   93  RCL 05   94  RCL 13   95  *   96  RCL 07 97  RCL 12   98  *   99  RCL 09 100  X^2 101  - 102  RCL 14 103  * 104  + 105  RCL 08 106  RCL 09 107  * 108  RCL 07 109  RCL 11         110  * 111  - 112  STO 16 113  RCL 15 114  * 115  + 116  RCL 00 117  / 118  STO 05 119  RCL 06 120  RCL 13 121  * 122  RCL 14 123  RCL 16 124  * 125  + 126  RCL 07 127  RCL 10 128  * 129  RCL 08 130  X^2 131  - 132  RCL 15 133  * 134  + 135  RCL 00 136  / 137  STO 06 138  RCL 05 139  RCL 04 140  RTN 141  LBL F 142  STO 00 143  XEQ IND 01 144  RCL 04 145  * 146  STO 16 147  RCL 00 148  XEQ IND 02 149  RCL 05 150  * 151  ST+ 16 152  RCL 00 153  XEQ IND 03 154  RCL 06 155  * 156  ST+ 16 157  RCL 16 158  RTN 159  GTO F 160  END

( 198 bytes / SIZE 017 )

Example:    Find 3 real numbers   a , b , c  such that  y = a + b x2 + c / ln x   best fits the following data:

 y 12.9141 0.4989 -18.2283 -40.5506 -68.2507 x 1.5 2 3 4 5

Then, estimate y(6) and y(7)

-Let's clear registers R07 thru R15:                        7.015  [CLRGX]  with an HP-41CX   ( or  [CLRG]  )
-We load these 3 functions into program memory:

01  LBL "T"         ( any global label, maximum of 6 characters )
02  1
03  RTN
04  LBL "U"        ( any global label, maximum of 6 characters )
05  X^2
06  RTN
07  LBL "V"        ( any global label, maximum of 6 characters )
08  LN
09  1/X
10  RTN

-We store their names into registers R01 R02 R03     "T" ASTO 01  "U"  ASTO 02  "V"  ASTO 03
-and we store the data ( GTO "LSF3" )

( enter y-values first )

12.9141   ENTER^    1.5   [A]      ( the Sigma+ key         x = 1.5  ends up in X-register
0.4989   ENTER^      2    [A]         in user mode )          x =   2   ends up in X-register
-18.2283   ENTER^     3    [A]
-40.5506   ENTER^     4    [A]                                                     ....  etc  ....
-68.2507   ENTER^     5    [A]

-Then  R/S  produces   2.4001               a =   2.4001 = R04
RDN               -3.0000               b = -3.0000 = R05
RDN                6.9999                c =  6.9999 = R06

Thus,    y = 2.4001 - 3.0000 x2 + 6.9999 / ln x

-To obtain y(6)       6   [F]  or   R/S   yields   -101.6933             6  is saved in R00
-To obtain y(7)       7   [F]  or   R/S   yields   -141.0029             7  is saved in R00

Note:

-No coefficient of determination is calculated by this program.

5°)  Linear Combination of n functions

Data Registers:               ( Registers    R00 R01 R02  ......  Rnn   are to be initialized and  R2n+1 thru Rn2+3n  are to be cleared )

•  R00 = n
•  R01 = f1 name      Rnn+1 = a1    R2n+1 = Sum f1.f1  .....................  Rn2+n+1 = Sum fn.f1      Rn2+2n+1  = Sum y.f1
•  R02 = f2 name      Rnn+2 = a2    R2n+2 = Sum f1.f2  .....................  Rn2+n+2  = Sum fn.f2     Rn2+2n+2  = Sum y.f2
............................................................................................................................................................................

•  Rnn = fn name      R2n = an        R3n   = Sum f1.fn  .....................     Rn2+2n  = Sum fn.fn       Rn2+3n   = Sum y.fn

Flags:  F01 is cleared by this program
Subroutines:   "LS"  synthetic version                  ( cf "Linear and non-linear systems for the HP-41" )
and the n functions  f1, f2, ... ,  fn  ( assuming x is in X-register upon entry )

-If you don't have an X-function module, replace lines 76 to 91 with:

ST+ X
E3
/
RCL 00
+
1
+
STO Z
SIGN
LBL 03
CLX
RCL IND Y
STO IND Z
ISG Y
CLX
ISG Z
GTO 03
LASTX
RTN

 01  LBL A   02  STO M   03  X<>Y   04  STO N   05  RCL 00   06  STO O   07  LBL 00   08  RCL IND O   09  RCL M   10  XEQ IND Y   11  RCL O   12  RCL 00   13  +   14  X<>Y   15  STO IND Y   16  LASTX   17  X^2   18  LASTX   19  +   20  ST+ Z   21  CLX   22  RCL N   23  * 24  ST+ IND Y   25  DSE O   26  GTO 00   27  RCL 00   28  X^2   29  LASTX   30  ST+ X   31  ST+ Y   32  RCL 00   33   E3   34  /   35  +   36  STO O   37  LBL 01   38  RCL IND X   39  RCL IND O   40  *   41  ST+ IND Z   42  RDN   43  DSE Y   44  DSE X   45  GTO 01   46  RCL 00 47  +   48  DSE O   49  GTO 01   50  RCL M   51  RTN   52  LBL "LSF"     53  3   54  RCL 00   55  ST+ Y   56  ST* Y   57   E2   58  /   59  +   60   E3   61  /   62  RCL 00   63  ST+ X   64  +   65  1   66  +   67  0   68  CF 01   69  XEQ "LS" 70  R^   71  STO 00   72  1   73  +   74  X^2   75  RCL 00    76  .1              77  %      78  +       79  1       80  +        81  STO Z      82   E3       83  /        84  +        85  REGMOVE   86  CLX      87  RCL 00     88   E3     89  /       90  +      91  RTN     92  LBL F 93  CLA     94  STO M     95  RCL 00   96  ST+ 00   97  STO N   98  LBL 02   99  RCL IND N 100  RCL M 101  XEQ IND Y 102  RCL IND 00 103  * 104  ST+ O 105  DSE 00 106  DSE N 107  GTO 02 108  X<> M 109  SIGN 110  X<> O 111  CLA 112  RTN 113  GTO F 114  END

( 178 bytes / SIZE n2+ 3n +1 )

Example:      Find 4 real numbers a1 , a2 , a3 , a4  such that    y = a1 + a2 sin 10.x + a3.x + a4 ln x   best fits the following data ( x in degrees ):

 y 5.22 6.69 7.28 7.01 3.06 -0.14 x 2 3 4 6 10 12

Then, evaluate y(5) and y(7)

-SIZE 029  or greater.
-Set your HP-41 in DEG mode.
-Clear registers R09 thru R28:                        9.028  [CLRGX]  with an HP-41CX   ( or  [CLRG]  )
-Load the 4 functions into program memory:

01  LBL "T"         ( any global label, maximum of 6 characters )
02  1
03  RTN
04  LBL "U"        ( any global label, maximum of 6 characters )
05  10
06  *
07  SIN
08  RTN
09  LBL "V"        ( any global label, maximum of 6 characters )
10  RTN
11  LBL "W"       ( any global label, maximum of 6 characters )
12  LN
13  RTN

-There are 4 functions:      4  STO 00
-Store their names into registers R01 thru R04     "T" ASTO 01    "U"  ASTO 02    "V"  ASTO 03   "W"  ASTO 04
-and store the data ( GTO "LSF" )

( enter y-values first )

5.22    ENTER^    2     [A]      ( the Sigma+ key       x = 2  ends up in X-register within  11 seconds
6.69    ENTER^    3     [A]         in user mode )        x = 3  ends up in X-register within  11 seconds
7.28    ENTER^    4     [A]
7.01    ENTER^    6     [A]                                                          ......  etc  ......
3.06    ENTER^   10    [A]
-0.14    ENTER^   12     [A]

-Then  R/S  or  XEQ "LSF"  produces ( in 31 seconds )        5.008  =  the control number of the solution:
a1 , a2 , a3 , a4  are in registers   R05 , R06 , R07 , R08
RCL 05  gives       a1    =  2.9957
RCL 06  -----       a2    =  4.0011
RCL 07  -----       a3    = -2.0014
RCL 08  -----       a4    =  7.0089

Thus,    y = 2.9957 + 4.0011 sin 10.x - 2.0014 x + 7.0089 ln x    is the best fitting curve for these data.

-To get  y(5)        5   [F]  or   R/S   yields   7.3339     ( in 3.5 seconds )      x = 5  is saved in L-register
-To get  y(7)        7   [F]  or   R/S   yields   6.3841     ( in 3.5 seconds )      x = 7  is saved in L-register

Remarks:

-Execution time is approximately proportional to n2
-No coefficient of correlation is calculated by this program.

6°)  Least-Squares Plane

-Given a set of n data points ( xi , yi , zi ) in space, we seek the best fitting plane  z = a x + b y + c

-We have to solve the 3x3 linear system

a Sum x2  + b Sum x.y + c Sum x  =  Sum x.z
a Sum x.y + b Sum y2  + c Sum y  =  Sum y.z
a Sum x    + b Sum y   + c n          =  Sum z

-The coefficient of determination is given by         r = [ a Sum x.z + b Sum y.z + c Sum z - (Sum z)2/n ] / [ Sum z2 - (Sum z)2/n ]

Data Registers:   R00 = Sum x     R02 = Sum y      R04 = Sum x.y    R06 = Sum z        R08 = Sum y.z         R10 = a    R11 = b    R12 = c
R01 = Sum x2    R03 = Sum y2    R05 = n               R07 = Sum x.z     R09 = Sum z2           R13 = r = correlation coefficient
Flags: /
Subroutines: /

 01  LBL I   02  SREG 04   03  CLS     04  SREG 00   05  CLS      06  RTN   07  LBL A   08  S+      09  X<> Z   10  ST+ 06   11  ST* Y   12  ST* L   13  ST* X   14  ST+ 09   15  X<>Y   16  ST+ 08   17  LASTX   18  ST+ 07   19  RCL 05   20  RTN   21  LBL "LR2"   22  RCL 01   23  RCL 05   24  *   25  RCL 00 26  X^2   27  -   28  STO 10   29  RCL 05   30  RCL 08   31  *   32  RCL 02   33  RCL 06   34  *   35  -   36  STO 12   37  *   38  STO 11   39  RCL 04   40  RCL 05   41  *   42  RCL 00       43  RCL 02   44  *   45  -   46  ST* 12   47  STO 13   48  RCL 05   49  RCL 07   50  * 51  RCL 00   52  RCL 06   53  *   54  -   55  *   56  ST- 11   57  LASTX   58  RCL 03   59  RCL 05   60  *   61  RCL 02   62  X^2   63  -   64  ST* 10   65  *   66  RCL 12       67  -   68  RCL 10   69  RCL 13   70  X^2   71  -   72  ST/ 11   73  /   74  STO 10   75  RCL 00 76  *   77  RCL 02   78  RCL 11   79  *   80  +   81  RCL 06   82  X<>Y   83  -   84  RCL 05   85  /   86  STO 12   87  RCL 06   88  *   89  RCL 07   90  RCL 10       91  *   92  +   93  RCL 08   94  RCL 11   95  *   96  +   97  RCL 09   98  RCL 06   99  X^2 100  RCL 05 101  / 102  ST- Z 103  - 104  / 105  STO 13 106  RCL 12 107  RCL 11 108  RCL 10 109  RTN 110  LBL F 111  RCL 10 112  * 113  X<>Y 114  RCL 11     115  * 116  + 117  RCL 12 118  + 119  RTN 120  GTO F 121  END

( 149 bytes / SIZE 014 )

 STACK INPUTS#i OUTPUTS#i OUTPUTS T / / r Z zi / c Y yi / b X xi i a

where   z = a x + b y + c  is the equation of the best fitting plane and  r = coefficient of correlation

Example:    Given the following data points

 zi 7 8 2 17 11 yi 2 1 0 3 2 xi 1 2 0 4 3

-First, we clear registers R00 thru R09   XEQ I  or  press  [I]  in user mode or  RTN  R/S
-We store the data

7   ENTER^        8   ENTER^       2   ENTER^      17  ENTER^      11  ENTER^
2   ENTER^        1   ENTER^       0   ENTER^       3   ENTER^       2   ENTER^
1   [A]                 2   [A]                0   [A]                4   [A]                3   [A]                       all in user mode

-Then press R/S  ( or XEQ "LR2" )  it yields

a =  2.425   = R10
RDN     b =  1.625   = R11
RDN     c =  1.55     = R12
RDN     r =  0.9822  = R13

So the least-squares plane is   z = 2.425 x + 1.625 y + 1.55

and the coefficient of determination is  0.9822  a relatively good fit!

-To obtain an estimation of, say  z(1;4) & z(7;3)   press

4  ENTER^   1  R/S  ( or [F]  in user mode )   z(1;4) =  10.475
3  ENTER^   7  R/S  ( or [F]  in user mode )    z(7;3) =  23.4
y  ENTER^   x  R/S  ( or [F]  in user mode )    ... and so on ...

7°)  Linear Combination of n functions ( Generalization to a Surface )

-Given a set of N data points ( xi , yi , zi ) in space, we seek a best fitting surface z = a1 f1(x,y) + a2 f2(x,y) + ...... + an fn(x,y)           ( n <= N )
where  f1 ,  f2 , ...... ,  fn      are n  known functions of 2 variables and  a1 , a2 , ...... , an   are unknown.
-The following program is quite similar to "LSF"

Data Registers:                    ( Registers    R00 R01 R02  ......  Rnn   are to be initialized and  R2n+1 thru Rn2+3n  are to be cleared )

•  R00 = n
•  R01 = f1 name      Rnn+1 = a1      R2n+1 = Sum f1.f1  .....................  Rn2+n+1 = Sum fn.f1      Rn2+2n+1  = Sum z.f1
•  R02 = f2 name      Rnn+2 = a2      R2n+2 = Sum f1.f2  .....................  Rn2+n+2  = Sum fn.f2     Rn2+2n+2  = Sum z.f2
............................................................................................................................................................................

•  Rnn = fn name       R2n = an          R3n   = Sum f1.fn  .....................     Rn2+2n  = Sum fn.fn      Rn2+3n   =  Sum z.fn

Flag:    F01 is cleared by this program
Subroutines:   "LS"  synthetic version             ( cf "Linear and non-linear systems for the HP-41" )
and the n functions  f1, f2, ... ,  fn  ( assuming x is in X-register and y is in Y-register upon entry )

-If you don't have an X-function module, replace lines 80 to 95 with:

ST+ X
E3
/
RCL 00
+
1
+
STO Z
SIGN
LBL 03
CLX
RCL IND Y
STO IND Z
ISG Y
CLX
ISG Z
GTO 03
LASTX
RTN

 01  LBL A   02  STO M   03  RDN   04  STO N   05  X<>Y   06  STO O   07  RCL 00   08  STO P   09  LBL 00   10  RCL IND P   11  RCL N   12  RCL M   13  XEQ IND Z      14  RCL P   15  RCL 00   16  +   17  X<>Y   18  STO IND Y   19  LASTX   20  X^2   21  LASTX   22  +   23  ST+ Z   24  CLX   25  RCL O 26  *   27  ST+ IND Y   28  DSE P   29  GTO 00   30  RCL 00   31  X^2   32  LASTX   33  ST+ X   34  ST+ Y   35  RCL 00   36   E3   37  /   38  +   39  STO O   40  LBL 01   41  RCL IND X      42  RCL IND O   43  *   44  ST+ IND Z   45  RDN   46  DSE Y   47  DSE X   48  GTO 01   49  RCL 00   50  + 51  DSE O   52  GTO 01   53  RCL N   54  RCL M   55  RTN   56  LBL "LSFXY"   57  3   58  RCL 00   59  ST+ Y   60  ST* Y   61   E2   62  /   63  +   64   E3   65  /   66  RCL 00   67  ST+ X   68  +   69  1   70  +   71  0   72  CF 01   73  XEQ "LS"   74  R^   75  STO 00 76  1   77  +   78  X^2   79  RCL 00     80  .1    81  %    82  +     83  1     84  +     85  STO Z    86   E3     87  /     88  +     89  REGMOVE      90  CLX      91  RCL 00    92   E3    93  /     94  +      95  RTN     96  LBL F   97  CLA    98  STO M   99  X<>Y 100  STO N 101  RCL 00 102  ST+ 00 103  STO P 104  LBL 02 105  RCL IND P 106  RCL N 107  RCL M 108  XEQ IND Z 109  RCL IND 00   110  * 111  ST+ O 112  DSE 00 113  DSE P 114  GTO 02 115  RCL N 116  RCL M 117  SIGN 118  X<> O 119  CLA 120  RTN 121  GTO F 122  END

( 195 bytes / SIZE n2+ 3n +1 )

Example:    Find 4 real numbers   a1 , a2 , a3 , a4  such that    y = a1 sin ( x + y ) + a2 ex / y + a3.x.y + a4 ln ( x.y )   best fits the following data ( x , y in radians )

 z -0.647 -6.301 -22.679 -57.46 -2.336 y 2 2 3 1 1 x 1 2 4 3 1

Then, evaluate z(1;4)

-SIZE 029  or greater.
-Clear registers R09 thru R28:                        9.028  [CLRGX]  with an HP-41CX   ( or  [CLRG]  )
-Load the 4 functions into program memory:

01  LBL "T"         ( any global label, maximum of 6 characters )
02  +
03  SIN
04  RTN
05  LBL "U"        ( any global label, maximum of 6 characters )
06  E^X
07  X<>Y
08  /
09  RTN
10  LBL "V"        ( any global label, maximum of 6 characters )
11  *
12  RTN
13  LBL "W"       ( any global label, maximum of 6 characters )
14  *
15  LN
16  RTN

-There are 4 functions:      4  STO 00
-Store their names into registers R01 thru R04     "T" ASTO 01    "U"  ASTO 02    "V"  ASTO 03    "W"   ASTO 04
-and store the data ( GTO "LSFXY" )

( enter z first , then y and x at last )

-0.647    ENTER^    2    ENTER^    1    [A]      ( the Sigma+ key      y = 2 and x = 1 end up in Y and X-registers
-6.301    ENTER^    2    ENTER^          [A]         in user mode )       y = 2 and x = 2 end up in Y and X-registers
-22.679   ENTER^    3    ENTER^    4    [A]
-57.460   ENTER^    1    ENTER^    3    [A]                                                  ......  etc  ......
-2.336    ENTER^    1    ENTER^          [A]

-Then  R/S  or  XEQ "LSFXY"    yields        5.008  =  the control number of the solution:         a1 , a2 , a3 , a4  are in registers   R05 , R06 , R07 , R08

RCL 05  gives       a1    =  2.0005
RCL 06  gives       a2    = -3.0000
RCL 07  gives       a3    =  3.9996
RCL 08  gives       a4    = -6.9985

Thus, the least-squares surface is    z = 2.0005 sin ( x + y ) - 3.0000 ex / y + 3.9996 x.y - 6.9985 ln ( x.y )

Then:   4  ENTER^ 1    R/S   ( or [F]  )   --------->    z(1;4) = 2.3393         ( y = 4 ends up in Y-register and x = 1 is saved in L-register )
... and so on ...

Remark:     If one function is, for instance,   sin [12.345 ( 2x + y )] , store 12.345 into an unused regiter ( say R99 ) and define this function by

ST+ X
+
RCL 99
*
SIN
RTN

Notes:

-No coefficient of determination  is calculated by this program.
-Status register P is used in LBL A and LBL F, therefore the fi subroutines must have no digit entry line
( unfortunately, register a wouldn't work because  XEQ IND Z  clears this register )

8°)  Least-Squares 3D-Hyperplane

-Given a set of n data points ( xi , yi , zi , ti ) in space, we seek the best fitting hyperplane  t = a x + b y + c z + d

-We have to solve the 4x4 linear system:

a Sum x2  + b Sum x.y + c Sum x.z + d Sum x = Sum x.t
a Sum x.y + b Sum y2  + c Sum y.z + d Sum y = Sum y.t
a Sum x.z + b Sum y.z + c Sum z2  + d Sum z =  Sum z.t
a Sum x    + b Sum y   + c Sum z   + d n         =  Sum t

-The correlation coefficient is given by         r = [ a Sum x.t + b Sum y.t + c Sum z.t + d Sum t - (Sum t)2/n ] / [ Sum t2 - (Sum t)2/n ]

Data Registers:   R00 = Sum x     R03 = Sum y2      R06 = Sum x.z   R09 = Sum y.t     R12 = Sum z2      R15 = a       R18 = d
R01 = Sum x2    R04 = Sum x.y     R07 = Sum x.t    R10 = Sum z.t     R13 = Sum t        R16 = b       R19 = r = coeff of correlation
R02 = Sum y     R05 = n                R08 = Sum y.z    R11 = Sum z       R14 = Sum t2      R17 = c       R20 thru R22: temp
Flags: /
Subroutines: /

-If you have a HP-41CX,  lines 02 to 07  may be replaced by   .014   CLRGX   SREG 00   RTN

 01  LBL I   02  SREG 09     03  CLS         04  SREG 06   05  CLS        06  SREG 00   07  CLS       08  RTN   09  LBL A   10  S+        11  LASTX   12  R^   13  ST+ 13   14  STO 15   15  ST* Y   16  RDN   17  ST+ 07   18  CLX   19  RCL 15   20  X<>Y   21  ST* Y   22  RDN   23  ST+ 09   24  RDN   25  ST+ 11   26  ST* 15   27  ST* L   28  ST* Z   29  ST* X   30  ST+ 12   31  X<> L   32  ST+ 06   33  RDN   34  X^2   35  ST+ 14   36  X<>Y   37  ST+ 08   38  RCL 15   39  ST+ 10   40  RCL 05   41  RTN   42  LBL "LR3"   43  RCL 01   44  RCL 05   45  * 46  RCL 00   47  X^2   48  -   49  RCL 05   50  RCL 06   51  *   52  RCL 00   53  RCL 11   54  *   55  -   56  STO 15   57  X^2   58  RCL 05   59  RCL 12   60  *   61  RCL 11   62  X^2   63  -   64  STO 16   65  ST/ 15   66  /   67  -   68  STO 19   69  RCL 04   70  RCL 05   71  *   72  RCL 00   73  RCL 02   74  *   75  -   76  RCL 05   77  RCL 08   78  *   79  RCL 02   80  RCL 11   81  *   82  -   83  STO 18   84  RCL 15   85  *   86  -   87  STO 20       88  RCL 05   89  RCL 07   90  * 91  RCL 00   92  RCL 13   93  *   94  -   95  RCL 05   96  RCL 10   97  *   98  RCL 11   99  RCL 13 100  * 101  - 102  STO 17 103  RCL 15 104  * 105  - 106  STO 21 107  RCL 03 108  RCL 05 109  * 110  RCL 02 111  X^2 112  - 113  RCL 18 114  X^2 115  RCL 16 116  / 117  - 118  RCL 05 119  RCL 09 120  * 121  RCL 02 122  RCL 13 123  * 124  - 125  RCL 17 126  RCL 18 127  * 128  RCL 16 129  ST/ 17 130  ST/ 18 131  / 132  - 133  STO 22     134  RCL 19 135  * 136  RCL 20 137  ST* 22 138  RCL 21 139  * 140  - 141  X<>Y 142  RCL 21 143  * 144  RCL 22 145  - 146  RCL 19 147  R^ 148  * 149  RCL 20 150  X^2 151  - 152  ST/ Z 153  / 154  ENTER^ 155  X<> 15 156  * 157  X<>Y 158  STO 16 159  RCL 18 160  * 161  + 162  ST- 17 163  RCL 13 164  RCL 00 165  RCL 15 166  * 167  - 168  RCL 02 169  RCL 16 170  * 171  - 172  RCL 11 173  RCL 17 174  * 175  - 176  RCL 05 177  / 178  STO 18     179  RCL 07 180  RCL 15 181  * 182  RCL 09 183  RCL 16 184  * 185  + 186  RCL 10 187  RCL 17 188  * 189  + 190  RCL 13 191  RCL 18 192  * 193  + 194  RCL 14 195  RCL 13 196  X^2 197  RCL 05 198  / 199  ST- Z 200  - 201  / 202  STO 19 203  SIGN 204  RCL 18 205  RCL 17 206  RCL 16 207  RCL 15 208  RTN 209  LBL F 210  RCL 15 211  * 212  X<>Y 213  RCL 16 214  * 215  + 216  X<>Y 217  RCL 17     218  * 219  + 220  RCL 18 221  + 222  RTN 223  GTO F 224  END

( 298 bytes / SIZE 023 )

 STACK INPUTS#i OUTPUTS#i OUTPUTS T ti / d Z zi / c Y yi / b X xi i a L / / r

where   t = a x + b y + c z + d   is the equation of the best-fitting hyperplane and  r = coefficient of correlation

Example:    Given the following data points

 ti 24 28 20 26 41 zi 1 2 -1 3 4 yi 2 3 4 0 1 xi 3 0 6 1 2

-First, we clear registers R00 thru R14     XEQ I  or  press  [I]  in user mode or  RTN  R/S
-We store the data

24  ENTER^       28  ENTER^      20  ENTER^      26  ENTER^      41  ENTER^
1   ENTER^        2   ENTER^      -1  ENTER^        3   ENTER^       4   ENTER^
2   ENTER^        3   ENTER^       4   ENTER^        0   ENTER^       1   ENTER^
3   [A]                 0   [A]                6   [A]                 1   [A]                2   [A]                       all in user mode

-Then press R/S  ( or XEQ "LR3" )  it yields

a =  2.1371   =   R15        ( in 5.6 seconds )
RDN     b =  3.8669   =  R16
RDN     c =  8.0766   =  R17
RDN     d =  0.3992   =  R18
LASTX   r =   0.9885

So the least-squares hyperplane is  t = 2.1371 x + 3.8669 y + 8.0766 z + 0.3992

and the correlation coefficient is   0.9885

-To obtain an estimation of     t(x,y,z)   press  z  ENTER^  y  ENTER^  x   R/S  ( or [F]  in user mode )
-For instance,

4  ENTER^
3  ENTER^
2    R/S       gives     t(2,3,4) = 48.5806

Note:

-There will be a DATA ERROR if all the z-values are equal
-In this case, change the name of the variables ( or use "LR2" ).

References:

[1]  F.Scheid - "Numerical Analysis" - Mc Graw Hill - ISBN  07-055197-9
[2]  Abramowitz and Stegun - "Handbook of Mathematical Functions" -  Dover Publications -  ISBN  0-486-61272-4
[3]  HP-41C  STAT PAC