hp41programs

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 ( x_i , y_i ) , we seek a best fitting curve y = a₁ f₁(x) + a₂ f₂(x) + ...... + a_n f_n(x)           ( n <= N )
  where  f₁ ,  f₂ , ...... ,  f_n      are n  known functions and  a₁ , a₂ , ...... , a_n   are unknown.
-The least-squares principle minimizes the sum of "vertical distances" to a curve:        SUM_i   [ y_i - a₁ f₁(x_i) + a₂ f₂(x_i) + ...... + a_n f_n(x_i) ]²
  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 x_i y_i - n µ_x µ_y ) / ( Sum x_i² - ( Sum x_i )² / n )                r = the linear correlation coefficient = ( Sum x_i y_i - n µ_x µ_y ) / (( n-1) s_x s_y)
     p =   µ_y - m.µ_x                                                                and       µ_x , µ_y ,  s_x , s_y    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 x² values               •   R03 = Summation of y² values                •   R05 = Number of data points

Flags: /
Subroutines: /
 
 

 
  ( 55 bytes / SIZE 006 )
 
 

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

  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 + A²/4 )^1/2             where  A = ( s_x/s_y - s_y/s_x ) / r    ;    r =  linear correlation coefficient  = ( Sum x_i y_i - n µ_x µ_y ) / (( n-1) s_x s_y)
     p =   µ_y - m.µ_x                                                                                    and       µ_x , µ_y ,  s_x , s_y    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 x² values               •   R03 = Summation of y² values                •   R05 = Number of data points

Flags: /
Subroutines: /
 
 

 
  ( 67 bytes / SIZE 006 )
 
 

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

  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 x² values           R03 = Summation of y² 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: /
 
 

 
  ( 109 bytes / SIZE 008 )
 
 

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

Example:
 

-Clear registers  R00 thru R05  ( for instance execute CLS )
-Then,  x_i  ENTER^   y_i  ENTER^    n_i   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 n_i's = 1 , lines 01 thru 22 are unuseful and the data may be stored by   y_i  ENTER^  x_i   S+   for all the points
 

2°)  Least-Squares Parabola
 

-Given a set of n data points ( x_i , y_i ) , we seek the best fitting parabola  y = a x² + b x + c

-We have to solve the 3x3 linear system

    a Sum x⁴ + b Sum x³ + c Sum x² = Sum x²y
    a Sum x³ + b Sum x² + c Sum x  = Sum x.y
    a Sum x² + b Sum x  + c n          = Sum y

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

Data Registers:   R00 = Sum x     R02 = Sum y      R04 = Sum x.y    R06 = Sum x³     R08 = Sum x²y        R12 = r
                              R01 = Sum x²    R03 = Sum y²    R05 = n               R07 = Sum x⁴     R09 = a , R10 = b , R11 = c
Flags: /
Subroutines: /
 
 

 
   ( 148 bytes / SIZE 013 )
 
 

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

Example:    Given the following data points
 

-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 x² + 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 ( x₀ , y₀ ) ,  ( x₁ , y₁ )  , ........................ ,  ( x_N , y_N ) ,
  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 x_i's are equally spaced
-We define orthogonal polynomials P_k,N(t)  by                                                                       ( t = (x-x₀)/h  with  h = x_i+1 - x_i )

   P_k,N(t) = SUM _i=0,1,...,k (-1)ⁱ  C_kⁱ C_k+iⁱ  t⁽ⁱ⁾/N⁽ⁱ⁾

     where   t⁽ⁱ⁾ = t(t-1)(t-2)....(t-i+1)  and  N⁽ⁱ⁾ = N(N-1)(N-2).....(N-i+1)      ( = 1 if  i = 0 )
       and    C_kⁱ = k!/(i!(k-i)!)

-The polynomial  p(x) can then be written  p(x) = a₀ P_0,N(t) + a₁ P_1,N(t) + .................... + a_m P_m,N(t)

    with   a_k = [ SUM _t=0,1,...,N  y_t P_k,N(t) ]/[ SUM _t=0,1,...,N  (P_k,N(t))² ]

-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:   c₀ + c₁ t + c₂ t² + ...... + c_m t^m

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

   P_-1,N(t) = 0 ,  P_0,N(t) = 1 ,  (k+1)(N-k) P_k+1,N(t) = (2k+1)(N-2t) P_k,N(t) - k(N+k+1) P_k-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 = y₀  •  R11 = y₁  •  R12 = y₂   .........................   •  R10+N = y_N

    >>>>    at the end          R11+N = a₀  ,  R12+N =  a₁  ,  .............................. ,  R11+N+m = a_m
                                     R12+N+m = c₀  , R13+N+m = c₁  , ........................... ,  R12+N+2m = c_m

                                     R13+N+2m to R13+N+3m = the coefficients of  P_m-1,N (t)
                                     R14+N+3m to R14+N+4m = the coefficients of  P_m,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
 
 

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

Example:
 

  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  c_i's            --- execution time = 1mn50s ---
                             RDN   17.021 = control number of the coefficients  a_i's

  and we get:

          p(t) = 5.771429 - 3.567857 P_1,N(t) - 1.005952 P_2,N(t) - 0.016667 P_3,N(t) + 0.037013 P_4,N(t)
          p(t) = 1.217965 + 2.088023 t + 0.093561 t² - 0.083586 t³ + 0.007197 t⁴

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 P_k,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 P_1,N(t) - 1.005952 P_2,N(t) - 0.016667 P_3,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 t² + 0.002778 t³   in 19 seconds

-However, the coefficient a₄ 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 ( x_i , y_i ), 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 )
 
 

 
   ( 198 bytes / SIZE 017 )
 

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

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 x² + 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 Rn²+3n  are to be cleared )

     •  R00 = n
     •  R01 = f₁ name      Rnn+1 = a₁    R2n+1 = Sum f₁.f₁  .....................  Rn²+n+1 = Sum f_n.f₁      Rn²+2n+1  = Sum y.f₁
     •  R02 = f₂ name      Rnn+2 = a₂    R2n+2 = Sum f₁.f₂  .....................  Rn²+n+2  = Sum f_n.f₂     Rn²+2n+2  = Sum y.f₂
         ............................................................................................................................................................................

     •  Rnn = f_n name      R2n = a_n        R3n   = Sum f₁.f_n  .....................     Rn²+2n  = Sum f_n.f_n       Rn²+3n   = Sum y.f_n

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  f₁, f₂, ... ,  f_n  ( 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
 
 

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

Example:      Find 4 real numbers a₁ , a₂ , a₃ , a₄  such that    y = a₁ + a₂ sin 10.x + a₃.x + a₄ ln x   best fits the following data ( x in degrees ):
 
 

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:
                                                                                                           a₁ , a₂ , a₃ , a₄  are in registers   R05 , R06 , R07 , R08
          RCL 05  gives       a₁    =  2.9957
          RCL 06  -----       a₂    =  4.0011
          RCL 07  -----       a₃    = -2.0014
          RCL 08  -----       a₄    =  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 n²
-No coefficient of correlation is calculated by this program.
 

6°)  Least-Squares Plane
 

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

-We have to solve the 3x3 linear system

    a Sum x²  + b Sum x.y + c Sum x  =  Sum x.z
    a Sum x.y + b Sum y²  + 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)²/n ] / [ Sum z² - (Sum z)²/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 x²    R03 = Sum y²    R05 = n               R07 = Sum x.z     R09 = Sum z²           R13 = r = correlation coefficient
Flags: /
Subroutines: /
 
 

 
  ( 149 bytes / SIZE 014 )
 
 

 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
 

-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 ( x_i , y_i , z_i ) in space, we seek a best fitting surface z = a₁ f₁(x,y) + a₂ f₂(x,y) + ...... + a_n f_n(x,y)           ( n <= N )
  where  f₁ ,  f₂ , ...... ,  f_n      are n  known functions of 2 variables and  a₁ , a₂ , ...... , a_n   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 Rn²+3n  are to be cleared )

     •  R00 = n
     •  R01 = f₁ name      Rnn+1 = a₁      R2n+1 = Sum f₁.f₁  .....................  Rn²+n+1 = Sum f_n.f₁      Rn²+2n+1  = Sum z.f₁
     •  R02 = f₂ name      Rnn+2 = a₂      R2n+2 = Sum f₁.f₂  .....................  Rn²+n+2  = Sum f_n.f₂     Rn²+2n+2  = Sum z.f₂
         ............................................................................................................................................................................

     •  Rnn = f_n name       R2n = a_n          R3n   = Sum f₁.f_n  .....................     Rn²+2n  = Sum f_n.f_n      Rn²+3n   =  Sum z.f_n

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  f₁, f₂, ... ,  f_n  ( 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
 
 

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

Example:    Find 4 real numbers   a₁ , a₂ , a₃ , a₄  such that    y = a₁ sin ( x + y ) + a₂ e^x / y + a₃.x.y + a₄ ln ( x.y )   best fits the following data ( x , y in radians )
 
 

Then, evaluate z(1;4)

-SIZE 029  or greater.
-XEQ RAD
-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:         a₁ , a₂ , a₃ , a₄  are in registers   R05 , R06 , R07 , R08

          RCL 05  gives       a₁    =  2.0005
          RCL 06  gives       a₂    = -3.0000
          RCL 07  gives       a₃    =  3.9996
          RCL 08  gives       a₄    = -6.9985

Thus, the least-squares surface is    z = 2.0005 sin ( x + y ) - 3.0000 e^x / 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 f_i 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 ( x_i , y_i , z_i , t_i ) 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 x²  + b Sum x.y + c Sum x.z + d Sum x = Sum x.t
    a Sum x.y + b Sum y²  + c Sum y.z + d Sum y = Sum y.t
    a Sum x.z + b Sum y.z + c Sum z²  + 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)²/n ] / [ Sum t² - (Sum t)²/n ]
 

Data Registers:   R00 = Sum x     R03 = Sum y²      R06 = Sum x.z   R09 = Sum y.t     R12 = Sum z²      R15 = a       R18 = d
                              R01 = Sum x²    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 t²      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
 
 

 
  ( 298 bytes / SIZE 023 )
 
 

 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

 
-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