Best Regression Fit
11-04-2017, 02:08 PM
Post: #1 Eddie W. Shore Senior Member Posts: 1,046 Joined: Dec 2013
Best Regression Fit
The program BESTFIT compares a set of regressions to determine a best fit. This simulates a feature presented on the Hewlett Packard HP 48S, HP 48G, HP 49G, and HP 50g. BESTFIT compares the correlations of the following four regression models:

1. Linear: y = a * x + b
2. Logarithmic: y = a * ln x + b
3. Exponential: y = b * a^x (y = b * e^(ln a * x))
4. Power: y = b * x^a

The output is a two element list: a string of the best fit equation and its corresponding correlation.

HP Prime Program: BESTFIT
Code:
 EXPORT BESTFIT(L1,L2) BEGIN // 2017-11-02 EWS // Simulate Best Fit // HP 48GX,49G,50g // initialize LOCAL clist,m,c,v,s,n; clist:={0,0,0,0}; // test // correlation is linear only // requires approx() // linear y=a*x+b clist:=approx(correlation(L1,L2)); // log y=a*LN(x)+b  c:=approx(correlation(LN(L1),L2)); IF IM(c)==0 THEN clist:=c; END; // exponential y=b*e^(a*x) c:=approx(correlation(L1,LN(L2))); IF IM(c)==0 THEN clist:=c; END; // power y=b*a^x c:=approx(correlation(LN(L1),LN(L2))); IF IM(c)==0 THEN clist:=c; END; // test // POS(L0^2,MAX(L0^2)) m:=POS(clist^2,MAX(clist^2)); c:=clist[m]; IF m==1 THEN v:=linear_regression(L1,L2); s:=STRING(v)+"+"+STRING(v)+ "*X"; RETURN {s,c}; END; IF m==2 THEN v:=logarithmic_regression(L1,L2); s:=STRING(v)+"+"+STRING(v)+ "*LN(X)"; RETURN {s,c}; END; IF m==3 THEN v:=exponential_regression(L1,L2); s:=STRING(v)+"*"+STRING(v)+ "^X"; RETURN {s,c}; END; IF m==4 THEN v:=power_regression(L1,L2); s:=STRING(v)+"*X^"+ STRING(v); RETURN {v,s}; END; END;

BESTFIT2: An Extended Version

Version 2 adds the following regressions:

5. Inverse: y = b + a/x
6. Simple Logistic: y = 1/(b + a*e^(-x))
7. Simple Quadratic: y = b + a*x^2
8. Square Root: y = √(a*x + b)

HP Prime Program: BESTFIT2
Code:
 EXPORT BESTFIT2(L1,L2) BEGIN // 2017-11-02 EWS // Simulate Best Fit // HP 48GX,49G,50g // additional models // initialize LOCAL clist,m,c,v,s; clist:={0,0,0,0,0,0,0,0}; // test // correlation is linear only // requires approx() // linear y=a*x+b clist:=approx(correlation(L1,L2)); // log y=a*LN(x)+b  c:=approx(correlation(LN(L1),L2)); IF IM(c)==0 THEN clist:=c; END; // exponential y=b*e^(a*x) c:=approx(correlation(L1,LN(L2))); IF IM(c)==0 THEN clist:=c; END; // power y=b*a^x c:=approx(correlation(LN(L1),LN(L2))); IF IM(c)==0 THEN clist:=c; END; // inverse y=b+a/x IF POS(L1,0)==0 THEN clist:=approx(correlation(1/L1, L2)); END; // simple logistic IF POS(L2,0)==0 THEN clist:=approx(correlation(e^(−L1), 1/L2)); END; // simple quadratic clist:=approx(correlation(L1^2, L2)); // square root  IF ΣLIST(L2≥0)==SIZE(L2) THEN clist:=approx(correlation(L1, L2^2)); END; // test // POS(L0^2,MAX(L0^2)) m:=POS(clist^2,MAX(clist^2)); c:=clist[m]; IF m==1 THEN v:=linear_regression(L1,L2); s:=STRING(v)+"+"+STRING(v)+ "*X"; END; IF m==2 THEN v:=logarithmic_regression(L1,L2); s:=STRING(v)+"+"+STRING(v)+ "*LN(X)"; END; IF m==3 THEN v:=exponential_regression(L1,L2); s:=STRING(v)+"*"+STRING(v)+ "^X"; END; IF m==4 THEN v:=power_regression(L1,L2); s:=STRING(v)+"*X^"+ STRING(v); END; // inverse IF m==5 THEN v:=linear_regression(1/L1,L2); s:=STRING(v)+"+"+STRING(v)+ "/X"; END; // simple logistic IF m==6 THEN v:=linear_regression(e^(−L1), 1/L2); s:="1/("+STRING(v)+"+"+ STRING(v)+"*e^(−X))"; END; // simple quadratic IF m==7 THEN v:=linear_regression(L1^2,L2); s:=STRING(v)+"+"+STRING(v)+ "*X^2"; END; // square root IF m==8 THEN v:=linear_regression(L1,L2^2); s:="√("+STRING(v)+"+"+ STRING(v)+"*X)"; END; RETURN {s,c}; END;

11-04-2017, 03:51 PM
Post: #2 salvomic Senior Member Posts: 1,366 Joined: Jan 2015
RE: Best Regression Fit
thanks Eddie,
good program.
It works better in Home: I tried in CAS with the second list with a series of LN() and it return "Error: Bad argument type", the same input is good in Home, however. Please, control it.

Salvo

∫aL√0mic (IT9CLU), HP Prime 50g 41CX 71b 42s 12C 15C - DM42 WP34s :: Prime Soft. Lib
12-17-2017, 04:17 PM
Post: #3
 akmon Member Posts: 188 Joined: Jun 2014
RE: Best Regression Fit
Very interesting program. I missed that feature on the HP Prime, till now.

I´ve done a little modification for better input data.

Code:
EXPORT BESTFIT2() BEGIN // 2017-11-02 EWS // Simulate Best Fit // HP 48GX,49G,50g // additional models // initialize EDITMAT(M1,{"INTRODUZCA LOS DATOS",{},{"X","Y"}}); L1:=mat2list(col(M1,1)); L2:=mat2list(col(M1,2)); LOCAL clist,m,c,v,s; clist:={0,0,0,0,0,0,0,0}; // test // correlation is linear only // requires approx() // linear y=a*x+b clist:=approx(correlation(L1,L2)); // log y=a*LN(x)+b  c:=approx(correlation(LN(L1),L2)); IF IM(c)==0 THEN clist:=c; END; // exponential y=b*e^(a*x) c:=approx(correlation(L1,LN(L2))); IF IM(c)==0 THEN clist:=c; END; // power y=b*a^x c:=approx(correlation(LN(L1),LN(L2))); IF IM(c)==0 THEN clist:=c; END; // inverse y=b+a/x IF POS(L1,0)==0 THEN clist:=approx(correlation(1/L1, L2)); END; // simple logistic IF POS(L2,0)==0 THEN clist:=approx(correlation(e^(−L1), 1/L2)); END; // simple quadratic clist:=approx(correlation(L1^2, L2)); // square root  IF ΣLIST(L2≥0)==SIZE(L2) THEN clist:=approx(correlation(L1, L2^2)); END; // test // POS(L0^2,MAX(L0^2)) m:=POS(clist^2,MAX(clist^2)); c:=clist[m]; IF m==1 THEN v:=linear_regression(L1,L2); s:=STRING(v)+"+"+STRING(v)+ "*X"; END; IF m==2 THEN v:=logarithmic_regression(L1,L2); s:=STRING(v)+"+"+STRING(v)+ "*LN(X)"; END; IF m==3 THEN v:=exponential_regression(L1,L2); s:=STRING(v)+"*"+STRING(v)+ "^X"; END; IF m==4 THEN v:=power_regression(L1,L2); s:=STRING(v)+"*X^"+ STRING(v); END; // inverse IF m==5 THEN v:=linear_regression(1/L1,L2); s:=STRING(v)+"+"+STRING(v)+ "/X"; END; // simple logistic IF m==6 THEN v:=linear_regression(e^(−L1), 1/L2); s:="1/("+STRING(v)+"+"+ STRING(v)+"*e^(−X))"; END; // simple quadratic IF m==7 THEN v:=linear_regression(L1^2,L2); s:=STRING(v)+"+"+STRING(v)+ "*X^2"; END; // square root IF m==8 THEN v:=linear_regression(L1,L2^2); s:="√("+STRING(v)+"+"+ STRING(v)+"*X)"; END; RETURN {s,c}; END;
12-27-2017, 03:31 PM (This post was last modified: 12-29-2017 03:52 AM by Namir.)
Post: #4 Namir Senior Member Posts: 690 Joined: Dec 2013
RE: Best Regression Fit
Here is my version of BESTFIT. It iterates over many powers for Y and X, including powers with fractions.

Code:
EXPORT BESTFIT(Mat, IdxX, IdxY) BEGIN // initialize L1:=mat2list(col(Mat, IdxX)); // X L2:=mat2list(col(Mat, IdxY)); // Y LOCAL Lx, Ly; LOCAL xpwr, ypwr, ix, iy, xscale, yscale; LOCAL xminIdx, xmaxIdx, yminIdx, ymaxIdx; LOCAL r2, coeffs; LOCAL bestR2, bestXpwr, bestYpwr; // combination of index ranges and scales should iterate // in range -4, -3.5, -3, ..., -0.5, 0, 0.5, 1, 1.5, ..., 4 // Adjust these values as you see fit xminIdx := -8; yminIdx := -8; xmaxIdx := 8; ymaxIdx := 8; xscale := 0.5; yscale := 0.5; bestR2 := -1; FOR iy FROM yminIdx TO ymaxIdx DO    // Transform Y values    ypwr := iy * yscale;     IF iy == 0 THEN      Ly := LN(L2);    ELSE       Ly := L2^ypwr;    END;            FOR ix FROM xminIdx TO xmaxIdx DO      xpwr := ix * xscale;         // Transform X values      IF ix == 0 THEN        Lx := LN(L1);      ELSE         Lx := L1^xpwr;      END;      r2 := approx(correlation(Lx, Ly))^2;      if r2 > bestR2 THEN        bestR2 := r2;        bestXpwr := xpwr;        bestYpwr := ypwr;        coeffs := linear_regression(Lx, Ly);        IF r2==1 THEN          RETURN {bestR2, bestYpwr, bestXpwr, coeffs, coeffs};        END;      END;   END; END; RETURN {bestR2, bestYpwr, bestXpwr, coeffs, coeffs}; END;

The function requires the name of a data matrix and the indices of the columns that select X and Y data. The function returns:

1) Best Rsqr value.
2) Best power for Y. Please note that a 0 means ln(x).
3) Best power for X. Please note that a 0 means ln(x).
5) Slope for best curve.
6) Intercept of best curve.

The fitted models are in the general form:

Y^ypwr = intercept + slope*X^xpwr

When ypwr or xpwr is 0, the term represents a natural logarithm.

You can alter the range for the loops, the scales for X and Y as you see fit. Just make sure you don't try to raise negative values to negative powers or powers with fractions. If you have date the best thing to do is to scale their values to be in the range of 1 to 2, using:

x(i) = 1 + (x(i) - min(x))/(max(x) - min(x))
y(i) = 1 + (y(i) - min(y))/(max(y) - min(y))

Such a scaling ensures that all transformations will not generate a run time errors.

Enjoy

Namir
12-27-2017, 09:14 PM (This post was last modified: 12-29-2017 03:53 AM by Namir.)
Post: #5 Namir Senior Member Posts: 690 Joined: Dec 2013
RE: Best Regression Fit version 2
The next version of the best fit, BESTFIT2, tests the power model between X and Y first, before testing all other transformations. If you have negative values in either X or Y the do not use this function:

Code:
EXPORT BESTFIT2(Mat, IdxX, IdxY) BEGIN // initialize L1:=mat2list(col(Mat, IdxX)); // X L2:=mat2list(col(Mat, IdxY)); // Y LOCAL Lx, Ly; LOCAL xpwr, ypwr, ix, iy, xscale, yscale; LOCAL xminIdx, xmaxIdx, yminIdx, ymaxIdx; LOCAL r2, coeffs; LOCAL bestR2, bestXpwr, bestYpwr; // combination of index ranges and scales should iterate // in range -4, -3.5, -3, ..., -0.5, 0, 0.5, 1, 1.5, ..., 4 // Adjust these values as you see fit xminIdx := -8; yminIdx := -8; xmaxIdx := 8; ymaxIdx := 8; xscale := 0.5; yscale := 0.5; // Test power model first bestR2 := approx(correlation(LN(L1), LN(L2)))^2; bestXpwr := 0; bestYpwr := 0; coeffs := linear_regression(LN(L1), LN(L2));         FOR iy FROM yminIdx TO ymaxIdx DO    // Transform Y values    ypwr := iy * yscale;     IF iy == 0 THEN      Ly := LN(L2);    ELSE       Ly := L2^ypwr;    END;            FOR ix FROM xminIdx TO xmaxIdx DO      xpwr := ix * xscale;         // Transform X values      IF ix == 0 THEN        Lx := LN(L1);      ELSE         Lx := L1^xpwr;      END;      r2 := approx(correlation(Lx, Ly))^2;      IF r2 > bestR2 THEN        bestR2 := r2;        bestXpwr := xpwr;        bestYpwr := ypwr;        coeffs := linear_regression(Lx, Ly);        IF r2==1 THEN          RETURN {bestR2, bestYpwr, bestXpwr, coeffs, coeffs};        END;      END;   END; END; RETURN {bestR2, bestYpwr, bestXpwr, coeffs, coeffs}; END;

I wrote this modified version because while testing Y=X^2, I was getting 1/Y^2 = 1/X^4 as the best model. While the answer is theoretically correct, it is not in the desired form.

Namir
12-28-2017, 12:54 PM (This post was last modified: 12-29-2017 03:53 AM by Namir.)
Post: #6 Namir Senior Member Posts: 690 Joined: Dec 2013
RE: Best Regression Fit
Here is a third version, BESTFIT3, that maps the X and Y data to values in the range of (1, 2), making all transformations possible:

The function returns:
1) Best Rsqr value.
2) Power of best Y transformation (0 means ln(y)).
3) Power of best X transformation (0 means ln(x)).
4) Best slope.
5) Best intercept.
6) Minimum X value.
7) Maximum X value.
8) Minimum Y value.
9) Maximum Y value.

Make sure you use the powers, slope, intercept, minima, and maxima in estimating Yhat values:

Yhat' = (slope * ((X-xmin)/(xmax-xmin)+1)^xpwr + intercept)^(1/ypwr)
Yhat = ymin + (Yhat'-1)*(ymax-ymin)

if xpwr=0 and ywpr!=0 use:

Yhat' = (slope * ln((X-xmin)/(xmax-xmin)+1) + intercept)^(1/ypwr)
Yhat = ymin + (Yhat'-1)*(ymax-ymin)

if xpwr!=0 and ywpr=0 use:

Yhat' = exp(slope * ((X-xmin)/(xmax-xmin)+1)^xpwr + intercept)
Yhat = ymin + (Yhat'-1)*(ymax-ymin)

if both xpwr and ypwr are zero, use:

Yhat' = exp(slope * ln((X-xmin)/(xmax-xmin)+1) + intercept)
Yhat = ymin + (Yhat'-1)*(ymax-ymin)

Use ln() functions when xpwr s 0 and use exp() when the ypwr is 0.

Code:
EXPORT BESTFIT3(Mat, IdxX, IdxY) BEGIN LOCAL Lx, Ly; LOCAL xpwr, ypwr, ix, iy, xscale, yscale; LOCAL xminIdx, xmaxIdx, yminIdx, ymaxIdx; LOCAL r2, coeffs; LOCAL bestR2, bestXpwr, bestYpwr; LOCAL xmin, xmax, ymin, ymax; // initialize L1:=mat2list(col(Mat, IdxX)); // X L2:=mat2list(col(Mat, IdxY)); // Y xmin:=MIN(L1); xmax:=MAX(L1); ymin:=MIN(L2); ymax:=MAX(L2); // Map data into (1,2) range L1:=L1-xmin; L1:=L1/(xmax-xmin)+1; L2:=L2-ymin; L2:=L2/(ymax-ymin)+1; // combination of index ranges and scales should iterate // in range -4, -3.5, -3, ..., -0.5, 0, 0.5, 1, 1.5, ..., 4 // Adjust these values as you see fit xminIdx := -8; yminIdx := -8; xmaxIdx := 8; ymaxIdx := 8; xscale := 0.5; yscale := 0.5; // Test power model first bestR2 := approx(correlation(LN(L1), LN(L2)))^2; bestXpwr := 0; bestYpwr := 0; coeffs := linear_regression(LN(L1), LN(L2)); IF bestR2==1 THEN   RETURN {bestR2, 0, 0, coeffs, coeffs, xmin, xmax, ymin, ymax}; END;        FOR iy FROM yminIdx TO ymaxIdx DO    // Transform Y values    ypwr := iy * yscale;     IF iy == 0 THEN      Ly := LN(L2);    ELSE       Ly := L2^ypwr;    END;            FOR ix FROM xminIdx TO xmaxIdx DO      xpwr := ix * xscale;         // Transform X values      IF ix == 0 THEN        Lx := LN(L1);      ELSE         Lx := L1^xpwr;      END;      r2 := approx(correlation(Lx, Ly))^2;      IF r2 > bestR2 THEN        bestR2 := r2;        bestXpwr := xpwr;        bestYpwr := ypwr;        coeffs := linear_regression(Lx, Ly);        IF r2==1 THEN          RETURN {bestR2, bestYpwr, bestXpwr, coeffs, coeffs, xmin, xmax, ymin, ymax};        END;      END;   END; END; RETURN {bestR2, bestYpwr, bestXpwr, coeffs, coeffs, xmin, xmax, ymin, ymax}; END;
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)