Multiple Linear Regrssion Using Neural Networks - Namir - 2026-03-01 17:50strong>
I asked Claude AI to write me a program to implement Shallow Feed-forward Neural Network Regression. Neural networks usually do not give you the coefficients of the regression model they use. Instead, you give these NNs the matrix X and column vector Y, and also a column of Xv values to predict the y values for. Neural networks can model with acceptable errors very complex nonlinear regression models.
Claude AI generated the function mlrNN() and the test function testNN. This code should bring HP Prime calculations to the AI and Neural Networks era!!
Here is the implementation for function mlrNN:
Code: // ============================================================
// mlrNN() -- Shallow Feedforward Neural Network Regression
// for the HP Prime Graphing Calculator (PPL language)
// ============================================================
//
// Architecture
// ------------
// Input layer : m neurons (one per predictor column)
// Hidden layer: h neurons (h = m + 2), sigmoid activation
// Output layer: 1 neuron, linear (identity) activation
//
// Weight matrices (bias absorbed as an extra column)
// W1 : h x (m+1) input -> hidden (col m+1 = bias)
// W2 : 1 x (h+1) hidden -> output (col h+1 = bias)
//
// Training
// * Column-wise z-score normalisation of x and y
// * Xavier / Glorot uniform weight initialisation
// * Batch gradient descent with back-propagation
// * Learning rate with slow exponential decay
// * Multiple random restarts; best model kept
//
// Parameters
// ----------
// x : data matrix (n rows x m columns)
// y : response column vector (n rows x 1 column)
// xv : prediction row vector (1 row x m columns)
//
// Returns
// -------
// { yHat, RMSE }
//
// yHat : predicted response for xv (a scalar).
// Computed by forward-propagating the SINGLE
// new observation xv through the trained network
// and denormalising the output.
//
// RMSE : root mean square error on the TRAINING set.
// Computed by forward-propagating EVERY row
// of x through the trained network, denormalising
// each fitted value, and comparing with the
// original y:
//
// yh_i = network(x_i) for i = 1 .. n
// RMSE = sqrt( SUM( (yh_i - y_i)^2 ) / n )
//
// IMPORTANT: yHat (prediction for the NEW point xv)
// is NOT the same as yh (fitted training
// values used only for the RMSE calculation).
//
// Usage Example
// -------------
// x := [[1,2],[3,4],[5,6],[7,8]];
// y := [[5],[11],[17],[23]];
// xv := [[4,5]];
// mlrNN(x, y, xv) --> { yHat_value, RMSE_value }
//
// ============================================================
EXPORT mlrNN(x, y, xv)
BEGIN
LOCAL n, m, h, lst;
LOCAL nTr, nEp, lr0, decay;
LOCAL i, j, ep, tr;
LOCAL s, v, lr, curE, bestE;
LOCAL xMu, xSd, yMu, ySd;
LOCAL xN, yN, xA;
LOCAL W1, W2, bestW1, bestW2;
LOCAL Z1, A1, AA, yh, er;
LOCAL gW2, gA, gZ, gW1;
LOCAL xvN, Z1v, A1v, AAv;
LOCAL yHat, RMSE;
// ---- dimensions ----
lst := SIZE(x);
n := lst(1);
m := lst(2);
// ---- hyperparameters ----
h := m + 2; // hidden neurons
nTr := 5; // random restarts
nEp := 800; // epochs per restart
lr0 := 0.5; // initial learning rate
decay := 0.998; // LR decay each epoch
// ========================================================
// NORMALISE x COLUMNS (z-score)
// xN(i,j) = (x(i,j) - mean_j) / sd_j
// ========================================================
xMu := MAKEMAT(0, 1, m);
xSd := MAKEMAT(0, 1, m);
FOR j FROM 1 TO m DO
s := 0;
FOR i FROM 1 TO n DO
s := s + x(i, j);
END;
xMu(1, j) := s / n;
s := 0;
FOR i FROM 1 TO n DO
s := s + (x(i, j) - xMu(1, j))^2;
END;
xSd(1, j) := √(s / n);
IF xSd(1, j) < 1E-15 THEN
xSd(1, j) := 1;
END;
END;
xN := MAKEMAT(0, n, m);
FOR i FROM 1 TO n DO
FOR j FROM 1 TO m DO
xN(i, j) := (x(i, j) - xMu(1, j))
/ xSd(1, j);
END;
END;
// ========================================================
// NORMALISE y (z-score)
// yN(i) = (y(i) - yMu) / ySd
// ========================================================
yMu := 0;
FOR i FROM 1 TO n DO
yMu := yMu + y(i, 1);
END;
yMu := yMu / n;
s := 0;
FOR i FROM 1 TO n DO
s := s + (y(i, 1) - yMu)^2;
END;
ySd := √(s / n);
IF ySd < 1E-15 THEN
ySd := 1;
END;
yN := MAKEMAT(0, n, 1);
FOR i FROM 1 TO n DO
yN(i, 1) := (y(i, 1) - yMu) / ySd;
END;
// ========================================================
// BUILD AUGMENTED INPUT xA = [xN | 1]
// xA is n x (m+1); last column = 1 (bias term)
// ========================================================
xA := MAKEMAT(0, n, m + 1);
FOR i FROM 1 TO n DO
FOR j FROM 1 TO m DO
xA(i, j) := xN(i, j);
END;
xA(i, m + 1) := 1;
END;
// ========================================================
// TRAINING WITH MULTIPLE RANDOM RESTARTS
// ========================================================
bestE := 1E30;
bestW1 := MAKEMAT(0, h, m + 1);
bestW2 := MAKEMAT(0, 1, h + 1);
FOR tr FROM 1 TO nTr DO
// ---- Xavier uniform initialisation ----
// W1: h x (m+1)
W1 := MAKEMAT(0, h, m + 1);
v := √(6 / (m + 1 + h));
FOR i FROM 1 TO h DO
FOR j FROM 1 TO m + 1 DO
W1(i, j) := (RANDOM - 0.5) * 2 * v;
END;
END;
// W2: 1 x (h+1)
W2 := MAKEMAT(0, 1, h + 1);
v := √(6 / (h + 1 + 1));
FOR j FROM 1 TO h + 1 DO
W2(1, j) := (RANDOM - 0.5) * 2 * v;
END;
lr := lr0;
// ---- Epoch loop ----
FOR ep FROM 1 TO nEp DO
// ==== FORWARD PASS ====
// Z1 = xA * W1^T (n x h)
Z1 := xA * TRN(W1);
// A1 = sigmoid(Z1) (n x h)
A1 := MAKEMAT(0, n, h);
FOR i FROM 1 TO n DO
FOR j FROM 1 TO h DO
s := Z1(i, j);
IF s > 15 THEN
A1(i, j) := 1;
ELSE
IF s < -15 THEN
A1(i, j) := 0;
ELSE
A1(i, j) := 1 / (1 + EXP(-s));
END;
END;
END;
END;
// AA = [A1 | 1] (n x (h+1))
AA := MAKEMAT(0, n, h + 1);
FOR i FROM 1 TO n DO
FOR j FROM 1 TO h DO
AA(i, j) := A1(i, j);
END;
AA(i, h + 1) := 1;
END;
// yh = AA * W2^T (n x 1)
yh := AA * TRN(W2);
// ==== TRAINING ERROR (on normalised data) ====
er := yh - yN; // (n x 1)
s := 0;
FOR i FROM 1 TO n DO
s := s + er(i, 1)^2;
END;
curE := √(s / n);
// Early stop if already very small
IF curE < 1E-12 THEN
BREAK;
END;
// ==== BACKPROPAGATION ====
// ---- Output-layer gradient ----
// gW2 = (1/n) * er^T * AA (1 x (h+1))
gW2 := TRN(er) * AA / n;
// ---- Hidden-layer gradient ----
// Error back-propagated to hidden activations
// gA(i,j) = er(i) * W2(1,j) j=1..h
gA := MAKEMAT(0, n, h);
FOR i FROM 1 TO n DO
FOR j FROM 1 TO h DO
gA(i, j) := er(i, 1) * W2(1, j);
END;
END;
// Through sigmoid derivative: a*(1-a)
// gZ(i,j) = gA(i,j) * A1(i,j) * (1 - A1(i,j))
gZ := MAKEMAT(0, n, h);
FOR i FROM 1 TO n DO
FOR j FROM 1 TO h DO
gZ(i, j) := gA(i, j)
* A1(i, j) * (1 - A1(i, j));
END;
END;
// gW1 = (1/n) * gZ^T * xA (h x (m+1))
gW1 := TRN(gZ) * xA / n;
// ==== UPDATE WEIGHTS ====
W1 := W1 - lr * gW1;
W2 := W2 - lr * gW2;
// Decay learning rate
lr := lr * decay;
END; // ep
// ---- Keep best model across restarts ----
IF curE < bestE THEN
bestE := curE;
bestW1 := W1;
bestW2 := W2;
END;
END; // tr
// ========================================================
// PREDICT yHat FOR THE NEW POINT xv
// (This is NOT the same as the training fitted values)
// ========================================================
// Normalise and augment xv -> xvN (1 x (m+1))
xvN := MAKEMAT(0, 1, m + 1);
FOR j FROM 1 TO m DO
xvN(1, j) := (xv(1, j) - xMu(1, j))
/ xSd(1, j);
END;
xvN(1, m + 1) := 1;
// Hidden layer
Z1v := xvN * TRN(bestW1); // 1 x h
A1v := MAKEMAT(0, 1, h);
FOR j FROM 1 TO h DO
s := Z1v(1, j);
IF s > 15 THEN
A1v(1, j) := 1;
ELSE
IF s < -15 THEN
A1v(1, j) := 0;
ELSE
A1v(1, j) := 1 / (1 + EXP(-s));
END;
END;
END;
// Augment with bias
AAv := MAKEMAT(0, 1, h + 1);
FOR j FROM 1 TO h DO
AAv(1, j) := A1v(1, j);
END;
AAv(1, h + 1) := 1;
// Output (normalised)
yHat := AAv * TRN(bestW2); // 1 x 1
// Denormalise
yHat := yHat(1, 1) * ySd + yMu;
// ========================================================
// COMPUTE RMSE ON TRAINING DATA
// Forward-propagate ALL rows of x through the best
// model, denormalise, and compare with original y.
// RMSE = sqrt( sum( (yh_i - y_i)^2 ) / n )
// ========================================================
// Forward pass using xA (already built)
Z1 := xA * TRN(bestW1); // n x h
A1 := MAKEMAT(0, n, h);
FOR i FROM 1 TO n DO
FOR j FROM 1 TO h DO
s := Z1(i, j);
IF s > 15 THEN
A1(i, j) := 1;
ELSE
IF s < -15 THEN
A1(i, j) := 0;
ELSE
A1(i, j) := 1 / (1 + EXP(-s));
END;
END;
END;
END;
AA := MAKEMAT(0, n, h + 1);
FOR i FROM 1 TO n DO
FOR j FROM 1 TO h DO
AA(i, j) := A1(i, j);
END;
AA(i, h + 1) := 1;
END;
yh := AA * TRN(bestW2); // n x 1
// Denormalise yh and compute RMSE vs original y
s := 0;
FOR i FROM 1 TO n DO
v := yh(i, 1) * ySd + yMu;
s := s + (v - y(i, 1))^2;
END;
RMSE := √(s / n);
// ========================================================
// RETURN { yHat, RMSE }
// ========================================================
RETURN {yHat, RMSE};
END;
And here is the code for the test program testNN:
Code: // ============================================================
// testNN() -- Test harness for the mlrNN() function
// for the HP Prime Graphing Calculator
// ============================================================
//
// Runs 6 test cases with known input-output relationships,
// calls mlrNN() on each, and reports yHat, RMSE, expected
// value, and absolute prediction error.
//
// Test Cases
// ----------
// 1. y = 2*x + 3 (simple linear, 1 predictor)
// 2. y = x^2 (quadratic, 1 predictor)
// 3. y = 3*x1 + 2*x2 + 1 (linear, 2 predictors)
// 4. y = x1*x2 (interaction, 2 predictors)
// 5. y = sin(x) (nonlinear periodic)
// 6. y = constant (7) (edge case)
//
// Each test prints:
// Expected value, yHat, prediction error, and RMSE.
// A PASS/FAIL verdict based on whether RMSE < threshold.
//
// Usage
// -----
// Just run testNN() from the command line or a program.
// Make sure mlrNN() is already stored on the calculator.
//
// ============================================================
EXPORT testNN()
BEGIN
LOCAL x, y, xv, result;
LOCAL yHat, RMSE, expect, errP;
LOCAL nPass, nFail, thr;
LOCAL i, j, k, v;
nPass := 0;
nFail := 0;
PRINT();
PRINT("================================");
PRINT(" mlrNN() Test Suite");
PRINT("================================");
PRINT("");
// ========================================================
// TEST 1: y = 2*x + 3 (simple linear)
// Training: x = 1..8, predict at x = 4.5
// Expected: 2*4.5 + 3 = 12
// ========================================================
PRINT("TEST 1: y = 2x + 3");
x := [[1],[2],[3],[4],[5],[6],[7],[8]];
y := MAKEMAT(0, 8, 1);
FOR i FROM 1 TO 8 DO
y(i, 1) := 2 * x(i, 1) + 3;
END;
xv := [[4.5]];
expect := 12;
thr := 2;
result := mlrNN(x, y, xv);
yHat := result(1);
RMSE := result(2);
errP := ABS(yHat - expect);
PRINT(" Expected : " + STRING(expect));
PRINT(" yHat : " + STRING(yHat, 2, 6));
PRINT(" PredErr : " + STRING(errP, 2, 6));
PRINT(" RMSE : " + STRING(RMSE, 2, 6));
IF RMSE < thr THEN
PRINT(" --> PASS");
nPass := nPass + 1;
ELSE
PRINT(" --> FAIL (RMSE >= "
+ STRING(thr) + ")");
nFail := nFail + 1;
END;
PRINT("");
// ========================================================
// TEST 2: y = x^2 (quadratic)
// Training: x = -5..5, predict at x = 3.5
// Expected: 12.25
// ========================================================
PRINT("TEST 2: y = x^2");
x := MAKEMAT(0, 11, 1);
y := MAKEMAT(0, 11, 1);
FOR i FROM 1 TO 11 DO
x(i, 1) := i - 6;
y(i, 1) := x(i, 1)^2;
END;
xv := [[3.5]];
expect := 12.25;
thr := 4;
result := mlrNN(x, y, xv);
yHat := result(1);
RMSE := result(2);
errP := ABS(yHat - expect);
PRINT(" Expected : " + STRING(expect));
PRINT(" yHat : " + STRING(yHat, 2, 6));
PRINT(" PredErr : " + STRING(errP, 2, 6));
PRINT(" RMSE : " + STRING(RMSE, 2, 6));
IF RMSE < thr THEN
PRINT(" --> PASS");
nPass := nPass + 1;
ELSE
PRINT(" --> FAIL (RMSE >= "
+ STRING(thr) + ")");
nFail := nFail + 1;
END;
PRINT("");
// ========================================================
// TEST 3: y = 3*x1 + 2*x2 + 1 (2-predictor linear)
// Training: 12 points, predict at (2.5, 3.5)
// Expected: 3*2.5 + 2*3.5 + 1 = 15.5
// ========================================================
PRINT("TEST 3: y = 3x1 + 2x2 + 1");
x := MAKEMAT(0, 12, 2);
y := MAKEMAT(0, 12, 1);
k := 0;
FOR i FROM 1 TO 4 DO
FOR j FROM 1 TO 3 DO
k := k + 1;
x(k, 1) := i;
x(k, 2) := j;
y(k, 1) := 3 * i + 2 * j + 1;
END;
END;
xv := [[2.5, 3.5]];
expect := 15.5;
thr := 2;
result := mlrNN(x, y, xv);
yHat := result(1);
RMSE := result(2);
errP := ABS(yHat - expect);
PRINT(" Expected : " + STRING(expect));
PRINT(" yHat : " + STRING(yHat, 2, 6));
PRINT(" PredErr : " + STRING(errP, 2, 6));
PRINT(" RMSE : " + STRING(RMSE, 2, 6));
IF RMSE < thr THEN
PRINT(" --> PASS");
nPass := nPass + 1;
ELSE
PRINT(" --> FAIL (RMSE >= "
+ STRING(thr) + ")");
nFail := nFail + 1;
END;
PRINT("");
// ========================================================
// TEST 4: y = x1 * x2 (interaction / product)
// Training: 4x4 grid on [1..4]x[1..4]
// Predict at (2.5, 3.5)
// Expected: 2.5 * 3.5 = 8.75
// ========================================================
PRINT("TEST 4: y = x1 * x2");
x := MAKEMAT(0, 16, 2);
y := MAKEMAT(0, 16, 1);
k := 0;
FOR i FROM 1 TO 4 DO
FOR j FROM 1 TO 4 DO
k := k + 1;
x(k, 1) := i;
x(k, 2) := j;
y(k, 1) := i * j;
END;
END;
xv := [[2.5, 3.5]];
expect := 8.75;
thr := 3;
result := mlrNN(x, y, xv);
yHat := result(1);
RMSE := result(2);
errP := ABS(yHat - expect);
PRINT(" Expected : " + STRING(expect));
PRINT(" yHat : " + STRING(yHat, 2, 6));
PRINT(" PredErr : " + STRING(errP, 2, 6));
PRINT(" RMSE : " + STRING(RMSE, 2, 6));
IF RMSE < thr THEN
PRINT(" --> PASS");
nPass := nPass + 1;
ELSE
PRINT(" --> FAIL (RMSE >= "
+ STRING(thr) + ")");
nFail := nFail + 1;
END;
PRINT("");
// ========================================================
// TEST 5: y = sin(x) (nonlinear periodic)
// Training: 13 points in [-3..3] step 0.5
// Predict at x = 1.0
// Expected: sin(1) ~ 0.8415
// ========================================================
PRINT("TEST 5: y = sin(x)");
x := MAKEMAT(0, 13, 1);
y := MAKEMAT(0, 13, 1);
FOR i FROM 1 TO 13 DO
v := (i - 7) * 0.5;
x(i, 1) := v;
y(i, 1) := SIN(v);
END;
xv := [[1.0]];
expect := SIN(1.0);
thr := 0.3;
result := mlrNN(x, y, xv);
yHat := result(1);
RMSE := result(2);
errP := ABS(yHat - expect);
PRINT(" Expected : " + STRING(expect, 2, 6));
PRINT(" yHat : " + STRING(yHat, 2, 6));
PRINT(" PredErr : " + STRING(errP, 2, 6));
PRINT(" RMSE : " + STRING(RMSE, 2, 6));
IF RMSE < thr THEN
PRINT(" --> PASS");
nPass := nPass + 1;
ELSE
PRINT(" --> FAIL (RMSE >= "
+ STRING(thr) + ")");
nFail := nFail + 1;
END;
PRINT("");
// ========================================================
// TEST 6: y = 7 (constant -- edge case)
// Training: 5 points with varying x, all y = 7
// Predict at x = 99
// Expected: 7
// ========================================================
PRINT("TEST 6: y = 7 (constant)");
x := [[1],[2],[3],[4],[5]];
y := [[7],[7],[7],[7],[7]];
xv := [[99]];
expect := 7;
thr := 1;
result := mlrNN(x, y, xv);
yHat := result(1);
RMSE := result(2);
errP := ABS(yHat - expect);
PRINT(" Expected : " + STRING(expect));
PRINT(" yHat : " + STRING(yHat, 2, 6));
PRINT(" PredErr : " + STRING(errP, 2, 6));
PRINT(" RMSE : " + STRING(RMSE, 2, 6));
IF RMSE < thr THEN
PRINT(" --> PASS");
nPass := nPass + 1;
ELSE
PRINT(" --> FAIL (RMSE >= "
+ STRING(thr) + ")");
nFail := nFail + 1;
END;
PRINT("");
// ========================================================
// SUMMARY
// ========================================================
PRINT("================================");
PRINT(" RESULTS: "
+ STRING(nPass) + " passed, "
+ STRING(nFail) + " failed (of 6)");
PRINT("================================");
RETURN {nPass, nFail};
END;
|