Least-squares line fitting: Difference between revisions

From JSXGraph Wiki
No edit summary
No edit summary
 
(3 intermediate revisions by the same user not shown)
Line 1: Line 1:
This little JXSGraph application finds the line or the circle which is the best fit for given set of points.
This little JXSGraph application finds the line or the circle which is the best fit for given set of points.
If the resulting line is green, it is a straight line. If the line is blue, it is a circle.


<jsxgraph width="600" height="600">
<jsxgraph width="600" height="600">
Line 37: Line 38:
brd.unsuspendUpdate();
brd.unsuspendUpdate();


bestFit(p1);
bestFit(p2);
var bestFit = function(p) {
//
//
// From here on, the best-fitting circle or line is found by least-squares fitting.
// bestFit, the best-fitting circle or line is found by least-squares fitting.
//
//
var i, j, r = [], rbar = [], x = [], y = [], z = [], A = [[0,0,0],[0,0,0],[0,0,0]], n, d,
var bestFit = function(p) {
    eigen, minIndex, minE, ev, c, xm, ym, zm, radius;
    var i, j, r = [], rbar = [], x = [], y = [], z = [], A = [[0,0,0],[0,0,0],[0,0,0]], n, d,
n = p.length;
        eigen, minIndex, minE, ev, c, xm, ym, zm, radius;
for (i=0;i<n;i++) {
    n = p.length;
    r.push([1.0, p[i].X(), p[i].Y()]);
    for (i=0;i<n;i++) {
    d = r[i][0]*r[i][0] + r[i][1]*r[i][1] + r[i][2]*r[i][2];
        r.push([1.0, p[i].X(), p[i].Y()]);
    r[i][0] = 1.0 - r[i][0]/d;
        d = r[i][0]*r[i][0] + r[i][1]*r[i][1] + r[i][2]*r[i][2];
    r[i][1] /= d;
        r[i][0] = 1.0 - r[i][0]/d;
    r[i][2] /= d;
        r[i][1] /= d;
}
        r[i][2] /= d;
    }


for (j=0;j<3;j++) {
    for (j=0;j<3;j++) {
    for (i=0,d=0;i<n;i++) {
        for (i=0,d=0;i<n;i++) {
         d += r[i][j];
            d += r[i][j];
        }
        d /= n;
         rbar[j] = d;
        for (i=0;i<n;i++) {
            r[i][j] -= d;
        }
     }
     }
    d /= n;
    rbar[j] = d;
     for (i=0;i<n;i++) {
     for (i=0;i<n;i++) {
         r[i][j] -= d;
         A[0][0] += r[i][0]*r[i][0];
        A[0][1] += r[i][0]*r[i][1];
        A[0][2] += r[i][0]*r[i][2];
        A[1][0] += r[i][1]*r[i][0];
        A[1][1] += r[i][1]*r[i][1];
        A[1][2] += r[i][1]*r[i][2];
        A[2][0] += r[i][2]*r[i][0];
        A[2][1] += r[i][2]*r[i][1];
        A[2][2] += r[i][2]*r[i][2];
     }
     }
}
for (i=0;i<n;i++) {
    A[0][0] += r[i][0]*r[i][0];
    A[0][1] += r[i][0]*r[i][1];
    A[0][2] += r[i][0]*r[i][2];
    A[1][0] += r[i][1]*r[i][0];
    A[1][1] += r[i][1]*r[i][1];
    A[1][2] += r[i][1]*r[i][2];
    A[2][0] += r[i][2]*r[i][0];
    A[2][1] += r[i][2]*r[i][1];
    A[2][2] += r[i][2]*r[i][2];
}


eigen = JXG.Math.Numerics.Jacobi(A);
    eigen = JXG.Math.Numerics.Jacobi(A);
minIndex = 0;
    minIndex = 0;
minE = eigen[0][0][0];
    minE = eigen[0][0][0];
for (j=1;j<3;j++) {
    for (j=1;j<3;j++) {
    if (eigen[0][j][j]<minE) {
        if (eigen[0][j][j]<minE) {
        minIndex = j;
            minIndex = j;
        minE = eigen[0][j][j];
            minE = eigen[0][j][j];
        }
     }
     }
}
    ev = [eigen[1][0][minIndex],eigen[1][1][minIndex],eigen[1][2][minIndex]];
ev = [eigen[1][0][minIndex],eigen[1][1][minIndex],eigen[1][2][minIndex]];
    c = -(rbar[0]*ev[0]+rbar[1]*ev[1]+rbar[2]*ev[2]);
c = -(rbar[0]*ev[0]+rbar[1]*ev[1]+rbar[2]*ev[2]);


xm = -ev[1];
    xm = -ev[1];
ym = -ev[2];
    ym = -ev[2];
zm = 2.0*(c+ev[0]);
    zm = 2.0*(c+ev[0]);
//console.log(c, c+ev[0]);
    //console.log(c, c+ev[0]);


// If c is close to zero, the best fittting object is a line.
    // If c is close to zero, the best fittting object is a line.
// The best threshold parameter has yet to be determined.
    // The best threshold parameter has yet to be determined.
// At the moment it is set to 0.01.
    // At the moment it is set to 0.01.
if (Math.abs(c)<0.01) {
    if (Math.abs(c)<0.01) {
    brd.create('line',[zm,xm,ym], {strokeColor:'green'});
        brd.create('line',[zm,xm,ym], {strokeColor:'green'});
      
     }  else {
}  else {
        radius = Math.sqrt((xm*xm+ym*ym-2*c*zm)/(zm*zm));
    var radius = Math.sqrt((xm*xm+ym*ym-2*c*zm)/(zm*zm));
        brd.create('circle',[[zm,xm,ym],radius]);
    brd.create('circle',[[zm,xm,ym],radius]);
    }
}
}; // end of bestFit()


};
bestFit(p1);
bestFit(p2);


</jsxgraph>
</jsxgraph>
Line 115: Line 114:


// Experiments with lines and circles:
// Experiments with lines and circles:
if (false) {
 
    // Plot random points on a line disturbed by a random factor
// 1) Plot random points on a line disturbed by a random factor
     var i, p = [], angle, xr, yr, delta = 0.1;
     var i, p1 = [], angle, xr, yr, delta = 0.1;


     // Random points are constructed which lie roughly on a line
     // Random points are constructed which lie roughly on a line
Line 127: Line 126:
         yr = 10*(Math.random()-0.5);
         yr = 10*(Math.random()-0.5);
         xr = 0.*yr+delta*(Math.random()-0.5);
         xr = 0.*yr+delta*(Math.random()-0.5);
         p.push(brd.create('point',[xr, yr], {withLabel:false}));
         p1.push(brd.create('point',[xr, yr], {withLabel:false}));
     }
     }
} else {
 
    // Plot random points on a circle disturbed by a random factor
// 2) Plot random points on a circle disturbed by a random factor
     var i, p = [], angle, co, si, delta = 0.2;
     var i, p2 = [], angle, co, si, delta = 0.2;
   
   
     // Random points are constructed which lie roughly on a circle
     // Random points are constructed which lie roughly on a circle
Line 142: Line 141:
         co = 4*Math.cos(angle)+delta*(Math.random()-0.5);
         co = 4*Math.cos(angle)+delta*(Math.random()-0.5);
         si = 4*Math.sin(angle)+delta*(Math.random()-0.5);
         si = 4*Math.sin(angle)+delta*(Math.random()-0.5);
         p.push(brd.create('point',[co+2, si-1], {withLabel:false}));
         p2.push(brd.create('point',[co+2, si-1], {withLabel:false}));
     }
     }
}
brd.unsuspendUpdate();
brd.unsuspendUpdate();


//
//
// From here on, the best-fitting circle or line is found by least-squares fitting.
// bestFit, the best-fitting circle or line is found by least-squares fitting.
//
//
var i, j, r = [], rbar = [], x = [], y = [], z = [], A = [[0,0,0],[0,0,0],[0,0,0]], n, d,
var bestFit = function(p) {
    eigen, minIndex, minE, ev, c, xm, ym, zm, radius;
    var i, j, r = [], rbar = [], x = [], y = [], z = [], A = [[0,0,0],[0,0,0],[0,0,0]], n, d,
n = p.length;
        eigen, minIndex, minE, ev, c, xm, ym, zm, radius;
for (i=0;i<n;i++) {
    n = p.length;
    r.push([1.0, p[i].X(), p[i].Y()]);
    for (i=0;i<n;i++) {
    d = r[i][0]*r[i][0] + r[i][1]*r[i][1] + r[i][2]*r[i][2];
        r.push([1.0, p[i].X(), p[i].Y()]);
    r[i][0] = 1.0 - r[i][0]/d;
        d = r[i][0]*r[i][0] + r[i][1]*r[i][1] + r[i][2]*r[i][2];
    r[i][1] /= d;
        r[i][0] = 1.0 - r[i][0]/d;
    r[i][2] /= d;
        r[i][1] /= d;
}
        r[i][2] /= d;
    }


for (j=0;j<3;j++) {
    for (j=0;j<3;j++) {
    for (i=0,d=0;i<n;i++) {
        for (i=0,d=0;i<n;i++) {
         d += r[i][j];
            d += r[i][j];
        }
        d /= n;
         rbar[j] = d;
        for (i=0;i<n;i++) {
            r[i][j] -= d;
        }
     }
     }
    d /= n;
    rbar[j] = d;
     for (i=0;i<n;i++) {
     for (i=0;i<n;i++) {
         r[i][j] -= d;
         A[0][0] += r[i][0]*r[i][0];
        A[0][1] += r[i][0]*r[i][1];
        A[0][2] += r[i][0]*r[i][2];
        A[1][0] += r[i][1]*r[i][0];
        A[1][1] += r[i][1]*r[i][1];
        A[1][2] += r[i][1]*r[i][2];
        A[2][0] += r[i][2]*r[i][0];
        A[2][1] += r[i][2]*r[i][1];
        A[2][2] += r[i][2]*r[i][2];
    }
 
    eigen = JXG.Math.Numerics.Jacobi(A);
    minIndex = 0;
    minE = eigen[0][0][0];
    for (j=1;j<3;j++) {
        if (eigen[0][j][j]<minE) {
            minIndex = j;
            minE = eigen[0][j][j];
        }
     }
     }
}
     ev = [eigen[1][0][minIndex],eigen[1][1][minIndex],eigen[1][2][minIndex]];
for (i=0;i<n;i++) {
     c = -(rbar[0]*ev[0]+rbar[1]*ev[1]+rbar[2]*ev[2]);
     A[0][0] += r[i][0]*r[i][0];
 
    A[0][1] += r[i][0]*r[i][1];
     xm = -ev[1];
    A[0][2] += r[i][0]*r[i][2];
     ym = -ev[2];
     A[1][0] += r[i][1]*r[i][0];
     zm = 2.0*(c+ev[0]);
    A[1][1] += r[i][1]*r[i][1];
     //console.log(c, c+ev[0]);
     A[1][2] += r[i][1]*r[i][2];
     A[2][0] += r[i][2]*r[i][0];
     A[2][1] += r[i][2]*r[i][1];
     A[2][2] += r[i][2]*r[i][2];
}


eigen = JXG.Math.Numerics.Jacobi(A);
    // If c is close to zero, the best fittting object is a line.
minIndex = 0;
    // The best threshold parameter has yet to be determined.
minE = eigen[0][0][0];
    // At the moment it is set to 0.01.
for (j=1;j<3;j++) {
    if (Math.abs(c)<0.01) {
    if (eigen[0][j][j]<minE) {
        brd.create('line',[zm,xm,ym], {strokeColor:'green'});
         minIndex = j;
    }  else {
         minE = eigen[0][j][j];
         radius = Math.sqrt((xm*xm+ym*ym-2*c*zm)/(zm*zm));
         brd.create('circle',[[zm,xm,ym],radius]);
     }
     }
}
}; // end of bestFit()
ev = [eigen[1][0][minIndex],eigen[1][1][minIndex],eigen[1][2][minIndex]];
c = -(rbar[0]*ev[0]+rbar[1]*ev[1]+rbar[2]*ev[2]);


xm = -ev[1];
bestFit(p1);
ym = -ev[2];
bestFit(p2);
zm = 2.0*(c+ev[0]);
//console.log(c, c+ev[0]);


// If c is close to zero, the best fittting object is a line.
// The best threshold parameter has yet to be determined.
// At the moment it is set to 0.01.
if (Math.abs(c)<0.01) {
    brd.create('line',[zm,xm,ym], {strokeColor:'green'});
   
}  else {
    var radius = Math.sqrt((xm*xm+ym*ym-2*c*zm)/(zm*zm));
    brd.create('circle',[[zm,xm,ym],radius]);
}
</source>
</source>


[[Category:Examples]]
[[Category:Examples]]
[[Category:Statistics]]
[[Category:Statistics]]

Latest revision as of 18:16, 9 November 2010

This little JXSGraph application finds the line or the circle which is the best fit for given set of points. If the resulting line is green, it is a straight line. If the line is blue, it is a circle.

The underlying JavaScript code

var brd = JXG.JSXGraph.initBoard('jxgbox',{boundingbox:[-5,5,5,-5], keepaspectratio:true, axis:true});
brd.suspendUpdate();

// Experiments with lines and circles:

// 1) Plot random points on a line disturbed by a random factor
    var i, p1 = [], angle, xr, yr, delta = 0.1;

    // Random points are constructed which lie roughly on a line
    // defined by y = 0.3*x+1.
    // delta*0.5 is the maximal distance in y-direction of the random
    // points from the line.
    brd.suspendUpdate();
    for (i=0;i<100;i++) {
        yr = 10*(Math.random()-0.5);
        xr = 0.*yr+delta*(Math.random()-0.5);
        p1.push(brd.create('point',[xr, yr], {withLabel:false}));
    }

// 2) Plot random points on a circle disturbed by a random factor
    var i, p2 = [], angle, co, si, delta = 0.2;
 
    // Random points are constructed which lie roughly on a circle
    // of radius 4 having the origin as center.
    // delta*0.5 is the maximal distance in x- and y- direction of the random
    // points from the circle line.
    for (i=0;i<100;i++) {
        angle = Math.random()*2*Math.PI;
 
        co = 4*Math.cos(angle)+delta*(Math.random()-0.5);
        si = 4*Math.sin(angle)+delta*(Math.random()-0.5);
        p2.push(brd.create('point',[co+2, si-1], {withLabel:false}));
    }
brd.unsuspendUpdate();

//
// bestFit, the best-fitting circle or line is found by least-squares fitting.
//
var bestFit = function(p) {
    var i, j, r = [], rbar = [], x = [], y = [], z = [], A = [[0,0,0],[0,0,0],[0,0,0]], n, d,
        eigen, minIndex, minE, ev, c, xm, ym, zm, radius;
    n = p.length;
    for (i=0;i<n;i++) {
        r.push([1.0, p[i].X(), p[i].Y()]);
        d = r[i][0]*r[i][0] + r[i][1]*r[i][1] + r[i][2]*r[i][2];
        r[i][0] = 1.0 - r[i][0]/d;
        r[i][1] /= d;
        r[i][2] /= d;
    }

    for (j=0;j<3;j++) {
        for (i=0,d=0;i<n;i++) {
            d += r[i][j];
        }
        d /= n;
        rbar[j] = d;
        for (i=0;i<n;i++) {
            r[i][j] -= d;
        }
    }
    for (i=0;i<n;i++) {
        A[0][0] += r[i][0]*r[i][0];
        A[0][1] += r[i][0]*r[i][1];
        A[0][2] += r[i][0]*r[i][2];
        A[1][0] += r[i][1]*r[i][0];
        A[1][1] += r[i][1]*r[i][1];
        A[1][2] += r[i][1]*r[i][2];
        A[2][0] += r[i][2]*r[i][0];
        A[2][1] += r[i][2]*r[i][1];
        A[2][2] += r[i][2]*r[i][2];
    }

    eigen = JXG.Math.Numerics.Jacobi(A);
    minIndex = 0;
    minE = eigen[0][0][0];
    for (j=1;j<3;j++) {
        if (eigen[0][j][j]<minE) {
            minIndex = j;
            minE = eigen[0][j][j];
        }
    }
    ev = [eigen[1][0][minIndex],eigen[1][1][minIndex],eigen[1][2][minIndex]];
    c = -(rbar[0]*ev[0]+rbar[1]*ev[1]+rbar[2]*ev[2]);

    xm = -ev[1];
    ym = -ev[2];
    zm = 2.0*(c+ev[0]);
    //console.log(c, c+ev[0]);

    // If c is close to zero, the best fittting object is a line.
    // The best threshold parameter has yet to be determined.
    // At the moment it is set to 0.01.
    if (Math.abs(c)<0.01) {
        brd.create('line',[zm,xm,ym], {strokeColor:'green'});
    }  else {
        radius = Math.sqrt((xm*xm+ym*ym-2*c*zm)/(zm*zm));
        brd.create('circle',[[zm,xm,ym],radius]);
    }
}; // end of bestFit()

bestFit(p1);
bestFit(p2);