Least-squares line fitting: Difference between revisions

From JSXGraph Wiki
No edit summary
No edit summary
 
(7 intermediate revisions by the same user not shown)
Line 1: Line 1:
This little JXSGraph application finds the line - described by homogeneous coordinates [a,b,c] - that minimizes
This little JXSGraph application finds the line or the circle which is the best fit for given set of points.
:<math> \sum_{i=1}^n (ax_i+by_i+cz_i)^2.</math>
If the resulting line is green, it is a straight line. If the line is blue, it is a circle.
 
Coming soon...


<jsxgraph width="600" height="600">
<jsxgraph width="600" height="600">
var brd = JXG.JSXGraph.initBoard('jxgbox',{boundingbox:[-5,5,5,-5], keepaspectratio:true, axis:true});
var brd = JXG.JSXGraph.initBoard('jxgbox',{boundingbox:[-5,5,5,-5], keepaspectratio:true, axis:true});
brd.suspendUpdate();
brd.suspendUpdate();
if (false) {
 
     var i, p = [], angle, xr, yr, delta = 0.1;
// 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
     // Random points are constructed which lie roughly on a line
Line 18: Line 19:
         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 {
 
     var i, p = [], angle, co, si, delta = 0.2;
// 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
     // Random points are constructed which lie roughly on a circle
Line 32: Line 34:
         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();


//
//
// Ab hier wird der beste Kreis, bzw. die beste Gerade ermittelt.
// bestFit, the best-fitting circle or line is found by least-squares fitting.
var 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;
var bestFit = function(p) {
n = p.length;
    var i, j, r = [], rbar = [], x = [], y = [], z = [], A = [[0,0,0],[0,0,0],[0,0,0]], n, d,
for (i=0;i<n;i++) {
        eigen, minIndex, minE, ev, c, xm, ym, zm, radius;
    r.push([1.0, p[i].X(), p[i].Y()]);
    n = p.length;
    d = r[i][0]*r[i][0] + r[i][1]*r[i][1] + r[i][2]*r[i][2];
    for (i=0;i<n;i++) {
    r[i][0] = 1.0 - r[i][0]/d;
        r.push([1.0, p[i].X(), p[i].Y()]);
    r[i][1] /= d;
        d = r[i][0]*r[i][0] + r[i][1]*r[i][1] + r[i][2]*r[i][2];
     r[i][2] /= d;
        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);
 
</jsxgraph>
 
===The underlying JavaScript code===
<source lang="javascript">
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}));
    }


for (j=0;j<3;j++) {
// 2) Plot random points on a circle disturbed by a random factor
     for (i=0,d=0;i<n;i++) {
    var i, p2 = [], angle, co, si, delta = 0.2;
         d += r[i][j];
    // 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}));
     }
     }
    d /= n;
brd.unsuspendUpdate();
     rbar[j] = d;
 
//
// 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++) {
     for (i=0;i<n;i++) {
         r[i][j] -= d;
         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 (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];
}
/*
for (i=0;i<3;i++) {
     for (j=0;j<3;j++) {
     for (j=0;j<3;j++) {
         //A[i][j] /= n;
         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);
    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 (Math.abs(c)<0.01) {
 
    brd.create('line',[zm,xm,ym], {strokeColor:'green'});
    // If c is close to zero, the best fittting object is a line.
      
    // The best threshold parameter has yet to be determined.
}  else {
    // At the moment it is set to 0.01.
    var radius = Math.sqrt((xm*xm+ym*ym-2*c*zm)/(zm*zm));
    if (Math.abs(c)<0.01) {
    brd.create('circle',[[zm,xm,ym],radius]);
        brd.create('line',[zm,xm,ym], {strokeColor:'green'});
}
     }  else {
</jsxgraph>
        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);


</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);