Best Regression Fit

[Namir]

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[1], coeffs[2], 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[1], coeffs[2], xmin, xmax, ymin, ymax};
       END;
     END;
  END;
END;
RETURN {bestR2, bestYpwr, bestXpwr, coeffs[1], coeffs[2], xmin, xmax, ymin, ymax};
END;