Statistical - Cubic Regression (Cubic Spline Fit) ?

[jonmoore]

When using desktop CAS/statistical packages, I often find that a cubic spline regression provides the best fit, for the data I'm analysing. The TI Nspire includes cubic regression as standard as do Casio calculators such as the CG50 (and its ancestors), and I wondered if there are any great examples of cubic regression on HP calculators.

For smaller data sets, I tend to use the Curve Fitting routines included in SandMath on my DM41x these days (originally found on the venerable AECROM), as this includes a best-fit routine but cubic regression\spline\curve isn't one of the options.

And any of the good routines I've encountered on RPL calculators 48/50 are emulations of the AECROM routines and the story is much the same on the Prime.

I'm sure something exists somewhere in the HP ecosystem. And it may be labelled polynomial regression as this is ultimately the expression describing the spline (be it Quadratic/Cubic/Quartic etc.). But I've not found it on my travels to date.

Any suggestions will be most gratefully received.

[Gjermund Skailand]

I suggest you take a look at hpcalc.
For the hp49/50
https://www.hpcalc.org/hp49/math/numeric/intpol23.zip
there is also a version for the hp48
br Gjermund

[jonmoore]

@Gjermund | Thanks for the recommendation. It's not exactly what I'm looking for, but it's something which will come in useful nonetheless.

I'm ideally looking for something that hooks into the statistical summation data set and then fits a cubic spline to the data points much as would happen with the existing regression options. The programs you link to can be used to achieve similar outcomes but there's a bit more work involved.

The thing that I'd ideally like to find is something that functions much like the old 41 AECROM curve fitting routines but with a cubic interpolation fit too. The Best Fit option of the AECROM routines is ideal as it samples the data against all the regression options and picks the best for the dataset (saving you the need to first output a scatter plot before manually picking a suitable regression option). It's unfortunate that its regression options doesn't include cubic, as the AECROM curve fitting routines rightly deserve their legendary status in the 41 community.

[rawi]

Hi Jonmore,

for the HP Prime you have the possibility to compute a cubic regression in the 2-variable statistics app. It is one of a total of 12 models you can select from to fit your data. See the manual for details.

[Massimo Gnerucci]

Did you try the Curve Fitting Module V2 by Ángel?
Among the rest it does include a "CSPLINE" (Cubic Spline Intpl.) program.

Available @ the usual site (TOS)

[jonmoore]

Splendid. The last two answers are exactly what I've been looking for.

Thanks @Massimo and @rawi.

And funnily enough, I already have Angel's Curve Fitting module on my DM41x, I simply forgot I'd put it there!

As to the Prime option, 'that'll learn me' for always turning to my 50g for non-CAS stuff such as statistical calculations.

[Thomas Klemm]

Derivation

With the Vandermonde matrix \(X\):

\(
X =
\begin{bmatrix}
1&x_{1}&x_{1}^{2}&x_{1}^{3} \\
1&x_{2}&x_{2}^{2}&x_{2}^{3} \\
1&x_{3}&x_{3}^{2}&x_{3}^{3} \\
\vdots &\vdots &\vdots &\vdots \\
1&x_{n}&x_{n}^{2}&x_{n}^{3} \\
\end{bmatrix}
\)

… and the corresponding \(y\) vector:

\(
y =
\begin{bmatrix}
y_{1} \\
y_{2} \\
y_{3} \\
\vdots\\
y_{n} \\
\end{bmatrix}
\)

… we can write the equation for the coefficients \(c\) as:

\(
X^\top X c = X^\top y
\)

This leads to:

\(
\begin{bmatrix}
n & \sum x & \sum x^2 & \sum x^3 \\
\sum x & \sum x^2 & \sum x^3 & \sum x^4 \\
\sum x^2 & \sum x^3 & \sum x^4 & \sum x^5 \\
\sum x^3 & \sum x^4 & \sum x^5 & \sum x^6 \\
\end{bmatrix}
\,
\begin{bmatrix}
c_{0} \\
c_{1} \\
c_{2} \\
c_{3} \\
\end{bmatrix}
=
\begin{bmatrix}
\sum y \\
\sum x y \\
\sum x^2 y \\
\sum x^3 y \\
\end{bmatrix}
\)

Since the matrix \(A = X^\top X \) is positive semidefinite, we can find a Cholesky decomposition:

\(
A = L \, L^\top
\)

Here the matrix \(L\) is a real lower triangular matrix with positive diagonal entries.

I dully followed the steps in the linked Lecture 7 Notes to calculate \(L\) and then used forward and back substitution to solve:

\(
L \, L^\top c = X^\top y
\)


Programs

CLΣ:
Clear stack and registers and set ΣREG to 20.

Σ+
Add a data point \((x, y)\).

Σ-
Remove a data point \((x, y)\).

CFIT
Calculate the coefficients of the cubic polynomial.

P(X)
Evaluate the cubic polynomial at \(x\).

This program is for the HP-42s, but can also be used with the HP-41C:

Code:

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

Registers

Matrix Registers:
00: \(t \to c_0\)
01: \(u \to c_1\)
02: \(v \to c_2\)
03: \(w \to c_3\)
04: \(n\)
05: \(a\)
06: \(b\)
07: \(c\)
08: \(b\)
09: \(c\)
10: \(d\)
11: \(d\)
12: \(e\)
13: \(f\)

Statistics Registers:
14: \(Σ x^3 \to c\)
15: \(Σ x^4 \to d\)
16: \(Σ x^5 \to e\)
17: \(Σ x^6 \to f\)
18: \(Σ x^2 y \to v\)
19: \(Σ x^3 y \to w\)
20: \(Σ x \to a\)
21: \(Σ x^2 \to b\)
22: \(Σ y \to t\)
23: \(Σ y^2\)
24: \(Σ x y \to u\)
25: \(n \to n\)

Correspondence to the matrix:
\(
\begin{bmatrix}
n & a & b & c & \, & t \\
& b & c & d & \, & u \\
& & d & e & \, & v \\
& & & f & \, & w \\
\end{bmatrix}

\begin{bmatrix}
04 & 05 & 06 & 07 & \, & 00 \\
& 08 & 09 & 10 & \, & 01 \\
& & 11 & 12 & \, & 02 \\
& & & 13 & \, & 03 \\
\end{bmatrix}
\)

Example

Let us find the cubic regression function for the following dataset:

(0, 1), (2, 0), (3, 3), (4, 5), (5, 4).

CLΣ
1 ENTER 0 Σ+
0 ENTER 2 Σ+
3 ENTER Σ+
5 ENTER 4 Σ+
4 ENTER 5 Σ+
CFIT

T: -0.3868
Z: 3.0687
Y: -5.0755
X: 0.9973

Therefore, the cubic regression function that best fits our data is:

\(
y = 0.9973 - 5.0755x + 3.0687x^2 - 0.3868x^3
\)

Let us estimate the value for \(x = 1\):

1 P(X)

-1.3962


References

• Polynomial Regression
• Vandermonde Matrix
• Proof for why a matrix multiplied by its transpose is positive semidefinite
• Cholesky Decomposition
• Positive Definite Systems
• Cubic Regression Calculator
  

[jonmoore]

Thanks for the math lesson, Thomas. The math wasn't the problem, it's simply that I like to be a lazy bugger at times, and if there's a unicorn solution available somewhere my lazy bones would rather achieve my goal with the least amount of key presses possible.

Having said all that, I'll definitely have a dabble with your program on both the 42 and 41. Always good to have options.

[Thomas Klemm]

Quote:achieve my goal with the least amount of key presses possible

You can do this by assigning the programs to dedicated keys on the HP-41C, or then creating a CUSTOM menu on the HP-42S.

Thank you for raising this question and for letting me learn something new today.
Now I might go a little deeper into the matter:

Quote:An alternative way to eliminate taking square roots in the \( \mathbf {LL} ^{\mathrm {*}} \) decomposition is to compute the Cholesky decomposition \( \mathbf {A} =\mathbf {LDL} ^{\mathrm {*}} \), then solving \( \mathbf {Ly} =\mathbf {b} \) for y, and finally solving \( \mathbf {DL} ^{\mathrm {*}}\mathbf {x} =\mathbf {y} \).

[jte]

I implemented cubic fit routines for the HP Prime’s Function app’s “sketch a function” feature, but those routines are only available via the “sketch a function” interface.

Ideally, numerical round-off should be taken into account when implementing these sorts of routines. As it was practical for this problem, I used an arithmetic with enough precision to allow for the calculations to be performed without any numerical round-off. (The input to these cubic fit routines are points with bounded integer coordinates — points from touch input. This limits the general applicability of the “sketch a function” cubic fit routines in the Prime.)

[rprosperi]

I have an HP User Library program for the HP-71B that I bought from the Library back in the 80's when I was enamored with cubic splines. If that's of interest, I can scan it and post a copy, just let me know. Same topic but sounds like it's not the same context that you're interested in.

[Werner]

I have the impression cubic regression (what Thomas described, and very well indeed) and cubic spline regression are mixed up.
Cubic spline regression typically takes 4 or 5 (depending on the number of sample points) splines, and determines the best fit for those splines (with of course endpoints and first and second derivatives equal). That is somewhat more involved than cubic regression, fitting a 3rd degree polynomial to statistical data, which would be the same as using a single spline.
Werner

[Gjermund Skailand]

This is a bit off topic as there is no regression involved,
https://people.maths.ox.ac.uk/trefethen/barycentric.pdf

but Barycentric Lagrange interpolation is a numerically stable interpolation method which fits to all points of (large) datasets and can be used for interpolation and also for derivatives.

https://www.hpcalc.org/hp49/math/numeric/chebhp50.zip
included in hpgcc3 for the hp50g only.
br Gjermund

[jonmoore]

(06-01-2022 02:31 AM)rprosperi Wrote:  I have an HP User Library program for the HP-71B that I bought from the Library back in the 80's when I was enamored with cubic splines. If that's of interest, I can scan it and post a copy, just let me know. Same topic but sounds like it's not the same context that you're interested in.

Thanks Bob. I don't have a HP-71B but I do have a Casio fx-880p. And it might be fun to attempt to convert the 71B program. Failing that, I have in the past managed to convert 71B programs to both my 28S and DM41x so it might make for a good way to while away a rainy Sunday afternoon.

I only recently got back into old school Basic handhelds and prices on 71B's are crazy high for us Brits these days - import duties and 20% VAT are the final nails in the coffin of that particular eBay transaction!

[jonmoore]

(06-01-2022 07:09 AM)Gjermund Skailand Wrote:  This is a bit off topic as there is no regression involved,
https://people.maths.ox.ac.uk/trefethen/barycentric.pdf

but Barycentric Lagrange interpolation is a numerically stable interpolation method which fits to all points of (large) datasets and can be used for interpolation and also for derivatives.

https://www.hpcalc.org/hp49/math/numeric/chebhp50.zip
included in hpgcc3 for the hp50g only.
br Gjermund

Thanks for the info Gjermund. This thread is turning into a great resource of the differing options.

[Thomas Klemm]

(06-01-2022 06:26 AM)Werner Wrote:  I have the impression cubic regression and cubic spline regression are mixed up.

We can use the same dataset as in the Cubic Regression Calculator: (0, 1), (2, 0), (3, 3), (4, 5), (5, 4)

Code:

0 1
2 0
3 3
4 5
5 4

Cubic spline interpolation

Equation

\(
f(x) = \begin{cases}3.1105 \cdot 10^{-1}\cdot x^3 -4.8828 \cdot 10^{-61}\cdot x^2 -1.7442\cdot x + 1.0000, & \text{if } x \in [0,2], \\-8.5465 \cdot 10^{-1}\cdot x^3 + 6.9942\cdot x^2 -1.5733 \cdot 10^{1}\cdot x + 1.0326 \cdot 10^{1}, & \text{if } x \in (2,3], \\-4.5930 \cdot 10^{-1}\cdot x^3 + 3.4360\cdot x^2 -5.0581\cdot x -3.4884 \cdot 10^{-1}, & \text{if } x \in (3,4], \\6.9186 \cdot 10^{-1}\cdot x^3 -1.0378 \cdot 10^{1}\cdot x^2 + 5.0198 \cdot 10^{1}\cdot x -7.4023 \cdot 10^{1}, & \text{if } x \in (4,5].\end{cases}
\)

We get piecewise defined cubic polynomials that fit all given points.

Evaluation

\(
f(1)=-0.43314
\)

Graph



Usage

Quote:Originally, spline was a term for elastic rulers that were bent to pass through a number of predefined points, or knots.
These were used to make technical drawings for shipbuilding and construction by hand, as illustrated in the figure.

Polynomial interpolation

Equation

\(
f(x)=1.6667 \cdot 10^{-2}\cdot x^{4}-5.6667 \cdot 10^{-1}\cdot x^{3}+ 3.6833\cdot x^{2}-5.7333\cdot x+ 1.0000
\)

We get a single polynomial that fits all given points.

Evaluation

\(
f(1)=-1.6000
\)

Graph



Usage

Quote:In numerical analysis, polynomial interpolation is the interpolation of a given data set by the polynomial of lowest possible degree that passes through the points of the dataset.

Conclusion

Use a cubic spline if you want a function that goes smoothly through the given data points.
Disadvantage: The function is defined piecewise.

Use polynomial interpolation if you want a single polynomial passing the given data points.
Disadvantage: Runge's phenomenon shows that for high values of n, the interpolation polynomial may oscillate wildly between the data points.

Use cubic regression when you want to get the cubic polynomial that best fits the given data points.
Disadvantage: The polynomial usually does not pass the given data points.

[jonmoore]

It really depends on the data sets you're working with. I often work with data that broadly speaking has an S-shaped characteristic (with outliers, but that needs a different treatment anyway) and cubic regression works very well with data of this nature.

But as ever with this type of scenario, YMMD.

Thanks for everybody's input so far.

[jonmoore]

As an aside, ref my use of the word spline.

I work a lot with 3d DCC packages (digital content creation packages such as Maya and Houdini) and CAD packages too (such as Rhino, which incidentally, started out life as a specialist tool for designing boat hulls), and splines are central to their workflows.

The beauty of working with cubic spines in particular is that they lead to wonderfully smooth curves. Bezier curves are another DCC curve type, but these are more typical of 2d creative packages such as Adobe's illustrator. They provide far more control but the nature of how they're defined make them less useful for statistical purposes.

The trick with splines in general is to use as few control points as possible. Using too many control points, especially if some control points are too closely spaced, can lead to unpredictable kinks. These kinks have a greater likelihood with cubic splines as they're not designed for rapid changes of direction.

Working with splines as design tools gives you a great intuition for how they function when applied to statistical data.

[Werner]

(05-31-2022 08:24 AM)jonmoore Wrote:  (..) I often find that a cubic spline regression provides the best fit, for the data I'm analysing.

You left out the one thing jonmoore was talking about, Thomas; cubic spline regression, fitting data to 3-4-5 cubic splines, depending on the number of data points (potentially hundreds, thousands).
Werner

[jonmoore]

(06-02-2022 06:01 AM)Werner Wrote:  

(05-31-2022 08:24 AM)jonmoore Wrote:  (..) I often find that a cubic spline regression provides the best fit, for the data I'm analysing.

You left out the one thing jonmoore was talking about, Thomas; cubic spline regression, fitting data to 3-4-5 cubic splines, depending on the number of data points (potentially hundreds, thousands).
Werner

A modern statistical analysis program delivers this functionality by choosing the appropriate control points from those available in the scatter plot so that the interpolated curve best describes the average of the data points.

Here's another reference source, in this case, a statistical analysis library for .NET.

https://www.extremeoptimization.com/Docu...lines.aspx