regression through zero
Is anyone aware of (or has already developed) source code that will fit a least-squares regression of specified order, AND can be forced through zero. I found the following program through Google that works really well but I'm unsure of how to modify it in a way that allows me to force the relationship through the origin. I don't understand the theory behind the gauss function (along with the author I guess)...
publicclass Regression{
// Apply least squares to raw data to determine the coefficients for
// an n-order equation: y = a0*X^0 + a1*X^1 + ... + an*X^n.
// Returns the coefficients for the solved equation, given a number
// of y and x data points. The rawData input is given in the form of
// {{y0, x0}, {y1, x1},...,{yn, xn}}.The coefficients returned by
// the regression are {a0, a1,...,an} which corresponds to
// {X^0, X^1,...,X^n}. The number of coefficients returned is the
// requested equation order (norder) plus 1.
staticdouble[] linear_equation(double rawData[][],int norder){
double a[][] =newdouble[norder+1][norder+1];
double b[] =newdouble[norder+1];
double term[] =newdouble[norder+1];
double ysquare = 0;
// step through each of the raw data entries
for (int i = 0; i < rawData.length; i++){
b[0] += rawData[i][0];//gives the sum of the y observations in b[0];
ysquare += rawData[i][0] * rawData[i][0];// sums the squares of the y-values.
// sum the x power values
double xpower = 1;//initializes xpower to value 1
for (int j = 0; j < norder+1; j++){
term[j] = xpower;
a[0][j] += xpower;
xpower = xpower * rawData[i][1];
}
// now set up the rest of rows in the matrix - multiplying each row by each term
for (int j = 1; j < norder+1; j++){
b[j] += rawData[i][0] * term[j];
for (int k = 0; k < b.length; k++){
a[j][k] += term[j] * term[k];
}
}
}
// solve for the coefficients
double coef[] = gauss(a, b);
// calculate the r-squared statistic
double ss = 0;
double yaverage = b[0] / rawData.length;
for (int i = 0; i < norder+1; i++){
double xaverage = a[0][i] / rawData.length;
ss += coef[i] * (b[i] - (rawData.length * xaverage * yaverage));
}
rsquared = ss / (ysquare - (rawData.length * yaverage * yaverage));
// solved the simultaneous equations via gauss
double [] newcoef =newdouble [coef.length + 1];
for (int i = 0; i < newcoef.length; i++){
System.arraycopy(coef,0,newcoef,0,i);
}
coef = newcoef;
coef[coef.length - 1] = rsquared;
return coef;
}
// it's been so long since I wrote this, that I don't recall the math
// logic behind it. IIRC, it's just a standard gaussian technique for
// solving simultaneous equations of the form: |A| = |B| * |C| where we
// know the values of |A| and |B|, and we are solving for the coefficients
// in |C|
staticdouble[] gauss(double ax[][],double bx[]){
double a[][] =newdouble[ax.length][ax[0].length];
double b[] =newdouble[bx.length];
double pivot;
double mult;
double top;
int n = b.length;
double coef[] =newdouble[n];
// copy over the array values - inplace solution changes values
for (int i = 0; i < ax.length; i++){
for (int j = 0; j < ax[i].length; j++){
a[i][j] = ax[i][j];
}
b[i] = bx[i];
}
for (int j = 0; j < (n-1); j++){
pivot = a[j][j];
for (int i = j+1; i < n; i++){
mult = a[i][j] / pivot;
for (int k = j+1; k < n; k++){
a[i][k] = a[i][k] - mult * a[j][k];
}
b[i] = b[i] - mult * b[j];
}
}
coef[n-1] = b[n-1] / a[n-1][n-1];
for (int i = n-2; i >= 0; i--){
top = b[i];
for (int k = i+1; k < n; k ++){
top = top - a[i][k] * coef[k];
}
coef[i] = top / a[i][i];
}
return coef;
}
// simple routine to print the equation for inspection
staticvoid print_equation(double coef[]){
for (int i = 0; i < coef.length; i++){
if (i == 0){
System.out.print("Y = ");
}else{
System.out.print(" + ");
}
System.out.print(coef[i] +"*X^" + i);
}
System.out.println("[r^2 = " + rsquared +"]");
}
}
Thanks for any assistance you can offer!
Andy
[9003 byte] By [
aliasa] at [2007-10-2 8:04:29]
Do you want linear regression through zero (i.e. just find a gradient) or (as with the source you posted) arbitrary polynomial regression?
I'll be using it to perform a series of simple regressions, and also a series of quadratic regressions. I actually don't need the simple regressions to be forced through zero, but the quadratics will definitely have to...Thanks,Andy
aliasa at 2007-7-16 21:57:49 >
Ideally I'd have a function that can achieve both linear and polynomial regression, AND required a boolean input to either force through the origin or not... I'm not averse to using two separate functions, one for the linear analyses, and one for the quadratic, but I haven't been able to find a reference to a polynomial function that includes the option of forcing through 0.
Thanks again.
aliasa at 2007-7-16 21:57:49 >
Compiled the code you posted, after adding the field it uses to pass data between linear_equation and print_equation. Seriously bad style. Not only that: I did a test with norder 2 and print_equation outputs a cubic! Have you really found this code helpful?
My trusty copy of CLR has given me sufficient insight that I should be able to post some working code to your spec today. I'm planning on using Jama to solve the matrix equation, because it's a lot better than any matrix solution code I can write given 2 days to do it.
/**
* Computes a least squares polynomial fit, optionally constraining the
* solution space to contain only polynomials which pass through the origin.
* @param data The data points, as <code>{{x_0, y_0}, {x_1, y_1},
* ...}</code>.
* @param polynomialOrder The order of the polynomial to fit.
* @param forceThroughZero Whether or not to constrain the polynomial to
* pass through the origin.
* @return A column matrix containing the coefficients of the polynomial. To
* extract them as a <code>double[]</code> use
* <code>getColumnPackedCopy()</code>.
*/
public static Matrix leastSquares(double[][] data,
int polynomialOrder,
boolean forceThroughZero)
{
// Let the basis functions f_i = x^i.
// Let the least square fit be F(x) = SUM_i=0^n c_i f_i(x) where
// n = polynomialOrder. If forceThroughZero then c_0 = 0, or
// alternatively n = polynomialOrder - 1 and f_i = x^(i + 1).
// Let matrices
//( f_0(x_0) ... f_n(x_0) )
//A = (.........)
//( f_0(x_m) ... f_n(x_m) )
// and
//( y_0 )
//y = ( ... )
//( y_m )
// where m = data.length - 1.
// Then the least squares constraint is equivalent (see CLR pp769f) to
//A^T A c = A^T y
int n = forceThroughZero ? polynomialOrder - 1 : polynomialOrder;
double[][] A_array = new double[data.length][n + 1];
double[][] y_array = new double[data.length][1];
for (int row = 0; row < A_array.length; row++)
{
double x_row = data[row][0];
A_array[row][0] = forceThroughZero ? x_row : 1;
for (int col = 1; col <= n; col++)
{
A_array[row][col] = A_array[row][col - 1] * x_row;
}
y_array[row][0] = data[row][1];
}
Matrix A = new Matrix(A_array);
Matrix y = new Matrix(y_array);
Matrix A_T = A.transpose();
Matrix A_Ty = A_T.times(y);
Matrix A_TA = A_T.times(A);
Matrix c = A_TA.solve(A_Ty);
// If !forceThroughZero c has the coefficients. If forceThroughZero we
// need to add in a 0 for the x^0 coefficient.
if (!forceThroughZero)
{
return c;
}
Matrix rv = new Matrix(polynomialOrder + 1, 1);
rv.setMatrix(1, polynomialOrder, 0, 0, c);
return rv;
}
Tested with Jama 1.0.2 ( http://math.nist.gov/javanumerics/jama/Jama-1.0.2.jar )
I'm extraordinarliy grateful for your help! I don't have time to test it right now (I'm a PhD in the physiological sciences, but I'll be sure to let you know how it works).Thanks!Andy
aliasa at 2007-7-16 21:57:49 >
It works really well, nothing to add.Cheers,Andy
aliasa at 2007-7-16 21:57:49 >
> I don't understand the theory behind the
> gauss function (along with the author I guess)...
>
> > // it's been so long since I wrote this, that I
> t I don't recall the math
> // logic behind it. IIRC, it's just a standard
> ard gaussian technique for
> // solving simultaneous equations of the form: |A|
> |A| = |B| * |C| where we
> // know the values of |A| and |B|, and we are
> are solving for the coefficients
>// in |C|
>
Now that the problem has been solved... that gauss function solves a linear system of equations, Ax=b. Its appears to be a numerically unstable implementation, so beware. With limited precision arithmetic Gaussian elimination requires "complete pivoting" or "partial pivoting" to reduce the error introduced when you pivot on a very small number.
Even with total pivoting it's not always the best way of solving matrix equations. That's why I use Jama: I've found in the past that it gives better solutions. I think it uses the properties of the matrix to select an appropriate decomposition.
> Even with total pivoting it's not always the best way
> of solving matrix equations. That's why I use Jama:
> I've found in the past that it gives better
> solutions. I think it uses the properties of the
> matrix to select an appropriate decomposition.
Hmm, interesting. I'll take a look at that. Thanks for the tip.
> > Even with total pivoting it's not always the best
> way
> > of solving matrix equations. That's why I use
> Jama:
> > I've found in the past that it gives better
> > solutions. I think it uses the properties of the
> > matrix to select an appropriate decomposition.
>
> Hmm, interesting. I'll take a look at that. Thanks
> for the tip.
I looked at some of the source for Jama, and it looks like it uses Gaussian elimination with partial pivoting to compute the LU decomposition. The other decomposition it uses is the QR decomposition which gives the least squares solution for overdetermined systems. As one would expect for a numerical library written in Java, the routines aren't optimized for any particular architecture. For this reason, a better matrix toolkit might be MTJ:
http://rs.cipr.uib.no/mtj/
which lets you plug in a native BLAS. I've found native BLAS's like ATLAS (http://math-atlas.sourceforge.net/)to be 4-8 times faster than my own naively written Java codes. Its a bit of a pain to get everything compiled and what not, but the performance difference is probably worth it if you do anything with large matrices.
Maybe my total pivoting was just numerically bad somehow. Or maybe decomposing into LU with pivoting is numerically better than straight elimination with pivoting.
Fantastic class, simple, elegant, easy to use. Works great for my test data, but when I throw large numbers at it I get decreasing accuracy.
First a medium size example with OK accuracy:
56,5609.982142857143,
112,2804.9910714285716,
168,1869.9940476190477,
224,1402.4955357142858,
280,1121.9964285714286,
336,934.9970238095239,
392,801.4260204081633,
448,701.2477678571429,
504,623.3313492063492
Plug this into the applet at http://www.wcrl.ars.usda.gov/cec/java2/polyrg1.htm
for an 8 factor polynomial and get a beautiful fit. When I use this class and then ask the resulting poly to predict the value at, say, point 280 I get
1121.9930701804597
pretty good results, but if we use larger numbers such as
56050,5604982.43794826,
112100,2802491.21897413,
168150,1868327.4793160867,
224200,1401245.609487065,
280250,1120996.4875896522,
336300,934163.7396580434,
392350,800711.7768497515,
448400,700622.8047435326,
504450,622775.8264386957
and ask it to predict the value at 336300 we get 961559.4288567305, which is starting to depart rather far from the desired answer. My hope was to use this with some cryptographically huge numbers, but that isn't going to work as-is.
What I was thinking was I might go into the Jama library and convert it from double to BigDecimal, but now I'm seeing it might be some subtle gausian thing that might be beyond my abilities... :(
If you're planning on using it for cryptographically large numbers, does that imply you're doing cryptography? If so, wouldn't you want to use a field other than doubles anyway?
My code uses BigDecimal, I haven't found a class that will perform a polynomial regression on BigDecimal so I was converting down to double while testing out the suitability of this class with smaller numbers. If the math works out then I'll invest the time to rework it for BigDecimal. Of course, at this point I don't know if the math isn't working out because of lost accuracy or if it's because of Gausian issues.
For those who want to get up and running with a quick test of this class, here is some simple code to do it. Here also is a simple method "polyResult" to find the result of whatever polynomial is returned.
public void testPolynomialMath() {
double[][] data = {{10.0, 1.0}, {20.0, 2.0}, {30.0, 4.0}, {40.0, 8.0}, {50.0, 16.0}};
final int FORTH_ORDER = 4;
double[] v = leastSquares(data, FORTH_ORDER).getColumnPackedCopy();
for (int i = 0; i < v.length; i++) {
System.out.println("test polynomial = " + v);
}
long testPoint = 40;
System.out.println("for test point " + testPoint + " the predicted answer is " + polyResult(testPoint,v));
}
private double polyResult(long x, double[] poly) {
double result = 0;
BigInteger x2 = new BigInteger(Long.toString(x));
for (int i = 0; i < poly.length; i++) {
result = result + poly * x2.pow(i).doubleValue();
}
return result;
}
A handy tip for evaluating polynomials: a_n x^n + ... + a_1 x + a_0 =
a_0 + x (a_1 + x (a_2 + ... (a_n-1 + x a_n) ...))So /**
* Evaluates a polynomial.
* @param coeff The coefficients a_i, beginning with the constant a_0.
* @param x The point at which to evaluate.
*/
public static double evalPoly(double[] coeff, double x)
{
double accum = 0;
for (int i = coeff.length - 1; i >= 0; i--)
{
accum = accum * x + coeff[i];
}
return accum;
}
If you want to adapt it to use BigInteger you can, but I don't see the point myself, unless you're expecting catastrophic loss of significance.
Hi, I really need help about the implementation of weighted regression. The above code is regression. Could anyone give some help in the method of weighted regression (each point with its own weight ). I do not know how to modify the above code to let it become weighted regression.
Any help from you will be very appreciated!!!
R.
I posted my weighted regression code here modified based on the original regression code. Please give me some suggestions if the modified regression code is correctly used to weighted regression. Or give me some testing data to let me test it. Thank for any help!!!
// Apply least squares to raw data to determine the coefficients for
// an n-order equation: y = a0*X^0 + a1*X^1 + ... + an*X^n.
// Returns the coefficients for the solved equation, given a number
// of y and x data points. The rawData input is given in the form of
// {{y0, x0}, {y1, x1},...,{yn, xn}}.The coefficients returned by
// the regression are {a0, a1,...,an} which corresponds to
// {X^0, X^1,...,X^n}. The number of coefficients returned is the
// requested equation order (norder) plus 1.
static double[] linear_equation(double rawData[][], int norder, double [] weights, int length) {
double a[][] = new double[norder+1][norder+1];
double b[] = new double[norder+1];
double term[] = new double[norder+1];
double ysquare = 0;
double weight = 0;
// step through each raw data entries
for (int i = 0; i < rawData.length; i++) {
// sum the y values
b[0] += rawData[i][0];
ysquare += rawData[i][0] * rawData[i][0];
weight += weights[i];
// sum the x power values
double xpower = 1;
for (int j = 0; j < norder+1; j++) {
term[j] = xpower;
a[0][j] += xpower;
xpower = xpower * rawData[i][1];
}
// now set up the rest of rows in the matrix - multiplying each row by each term
for (int j = 1; j < norder+1; j++) {
b[j] += rawData[i][0] * term[j];
for (int k = 0; k < b.length; k++) {
a[j][k] += term[j] * term[k];
}
}
}
// solve for the coefficients
double coef[] = gauss(a, b);
//Squares
double ss = 0;
double yaverage = b[0] / length;
double waverage = weight / length;
for (int i = 0; i < norder+1; i++) {
double xaverage = a[0][i] / length;
ss += coef[i] * (b[i] - (length * xaverage * yaverage));
}
double rsquared = ss / (waverage * (ysquare - (length * yaverage * yaverage)));
double [] newcoef = new double [coef.length + 1];
for (int i = 0; i < newcoef.length; i++){
System.arraycopy(coef,0,newcoef,0,i);
}
coef = newcoef;
coef[coef.length - 1] = rsquared;
return coef;
}
static double[] gauss(double ax[][], double bx[]) {
double a[][] = new double[ax.length][ax[0].length];
double b[] = new double[bx.length];
double pivot;
double mult;
double top;
int n = b.length;
double coef[] = new double[n];
// copy over the array values - inplace solution changes values
for (int i = 0; i < ax.length; i++) {
for (int j = 0; j < ax[i].length; j++) {
a[i][j] = ax[i][j];
}
b[i] = bx[i];
}
for (int j = 0; j < (n-1); j++) {
pivot = a[j][j];
for (int i = j+1; i < n; i++) {
mult = a[i][j] / pivot;
for (int k = j+1; k < n; k++) {
a[i][k] = a[i][k] - mult * a[j][k];
}
b[i] = b[i] - mult * b[j];
}
}
coef[n-1] = b[n-1] / a[n-1][n-1];
for (int i = n-2; i >= 0; i--) {
top = b[i];
for (int k = i+1; k < n; k ++) {
top = top - a[i][k] * coef[k];
}
coef[i] = top / a[i][i];
}
return coef;
}
hello allPlease could some one take a look at this thread http://forum.java.sun.com/thread.jspa?threadID=5171483&tstart=0I am relatively new to all these topicsRegardsPac
> hello all
>
> Please could some one take a look at this thread
>
> http://forum.java.sun.com/thread.jspa?threadID=5171483
> &tstart=0
>
> I am relatively new to all these topics
>
> Regards
> Pac
Jeez Pac, you cross posted your message three times and now you're also hijacking other threads?
I think I'll pass, thank you.
