Best Regression Fit

[Namir]

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