Best Regression Fit

[Namir]

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