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 ( 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: /

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 + 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: /

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 )

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:

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: /

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:

x_i	2	5	10	16	21
y_i	4	7	12	19	24
n_i	3	4	7	5	2

-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: /

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	y_i	/	b
X	x_i	i	a

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

Example: Given the following data points

y_i	0	3	4	5	1
x_i	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 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

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(a_i)
X	m	bbb.eee(c_i)

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 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 , youhave 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 )

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 x² + 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 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

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 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 ):

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:
                                                                                                           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: /

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	z_i	/	c
Y	y_i	/	b
X	x_i	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

z_i	7	8	2	17	11
y_i	2	1	0	3	2
x_i	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 ( 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

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 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 )

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

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

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	t_i	/	d
Z	z_i	/	c
Y	y_i	/	b
X	x_i	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

t_i	24	28	20	26	41
z_i	1	2	-1	3	4
y_i	2	3	4	0	1
x_i	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