js实现曲线拟合(最小二乘法)

写在前面

很多时候我们需要把一组点用平滑曲线连接,从点到曲线,就是曲线拟合。而最小二乘法是实现曲线拟合最简单的方法,虽然效果不是很好,但胜在简单

一.曲线拟合

用连续曲线近似地刻画或比拟平面上离散点组所表示的坐标之间的函数关系的一种数据处理方法

嗯,说白了就是用平滑曲线连接。手握铅笔很容易搞定,但编程实现有一定难度,我们下面使用的方法叫最小二乘法,至于原理、数学公式、证明过程,抱歉,我不知道。

二.完整代码

P.S.下面的代码是根据某Java实现改写而来的,点我查看原文

function PolyFitForcast() {
    /**
     * <p>
     * 函数功能:最小二乘法曲线拟合
     * </p>
     * <p>
     * 方程:Y = a(0) + a(1) * (X - X1)+ a(2) * (X - X1)^2 + ..... .+ a(m) * (X -
     * X1)^m X1为X的平均数
     * </p>
     * 
     * @param x
     *            实型一维数组,长度为 n. 存放给定 n 个数据点的 X 坐标
     * @param y
     *            实型一维数组,长度为 n.存放给定 n 个数据点的 Y 坐标
     * @param n
     *            变量。给定数据点的个数
     * @param a
     *            实型一维数组,长度为 m.返回 m-1 次拟合多项式的 m 个系数
     * @param m
     *            拟合多项式的项数,即拟合多项式的最高次数为 m-1. 要求 m<=n 且m<=20。若 m>n 或 m>20
     *            ,则本函数自动按 m=min{n,20} 处理.
     *            <p>
     *            Date:2007-12-25 16:21 PM
     *            </p>
     * @author qingbao-gao
     * @return 多项式系数存储数组
     */
    function PolyFit(x, y, n, a, m) {
        var i, j, k;
        var z, p, c, g, q = 0, d1, d2;
        var s = new Array(20);
        var t = new Array(20);
        var b = new Array(20);
        var dt = new Array(3);
        for (i = 0; i <= m - 1; i++) {
            a[i] = 0.0;
        }
        if (m > n) {
            m = n;
        }
        if (m > 20) {
            m = 20;
        }
        z = 0.0;
        for (i = 0; i <= n - 1; i++) {
            z = z + x[i] / (1.0 * n);
        }
        b[0] = 1.0;
        d1 = 1.0 * n;
        p = 0.0;
        c = 0.0;
        for (i = 0; i <= n - 1; i++) {
            p = p + (x[i] - z);
            c = c + y[i];
        }
        c = c / d1;
        p = p / d1;
        a[0] = c * b[0];
        if (m > 1) {
            t[1] = 1.0;
            t[0] = -p;
            d2 = 0.0;
            c = 0.0;
            g = 0.0;
            for (i = 0; i <= n - 1; i++) {
                q = x[i] - z - p;
                d2 = d2 + q * q;
                c = c + y[i] * q;
                g = g + (x[i] - z) * q * q;
            }
            c = c / d2;
            p = g / d2;
            q = d2 / d1;
            d1 = d2;
            a[1] = c * t[1];
            a[0] = c * t[0] + a[0];
        }
        for (j = 2; j <= m - 1; j++) {
            s[j] = t[j - 1];
            s[j - 1] = -p * t[j - 1] + t[j - 2];
            if (j >= 3)
                for (k = j - 2; k >= 1; k--) {
                    s[k] = -p * t[k] + t[k - 1] - q * b[k];
                }
            s[0] = -p * t[0] - q * b[0];
            d2 = 0.0;
            c = 0.0;
            g = 0.0;
            for (i = 0; i <= n - 1; i++) {
                q = s[j];
                for (k = j - 1; k >= 0; k--) {
                    q = q * (x[i] - z) + s[k];
                }
                d2 = d2 + q * q;
                c = c + y[i] * q;
                g = g + (x[i] - z) * q * q;
            }
            c = c / d2;
            p = g / d2;
            q = d2 / d1;
            d1 = d2;
            a[j] = c * s[j];
            t[j] = s[j];
            for (k = j - 1; k >= 0; k--) {
                a[k] = c * s[k] + a[k];
                b[k] = t[k];
                t[k] = s[k];
            }
        }
        dt[0] = 0.0;
        dt[1] = 0.0;
        dt[2] = 0.0;
        for (i = 0; i <= n - 1; i++) {
            q = a[m - 1];
            for (k = m - 2; k >= 0; k--) {
                q = a[k] + q * (x[i] - z);
            }
            p = q - y[i];
            if (Math.abs(p) > dt[2]) {
                dt[2] = Math.abs(p);
            }
            dt[0] = dt[0] + p * p;
            dt[1] = dt[1] + Math.abs(p);
        }
        return a;
    }// end

    /**
     * <p>
     * 对X轴数据节点球平均值
     * </p>
     * 
     * @param x
     *            存储X轴节点的数组
     *            <p>
     *            Date:2007-12-25 20:21 PM
     *            </p>
     * @author qingbao-gao
     * @return 平均值
     */
    function average(x) {
        var ave = 0;
        var sum = 0;
        if (x !== null) {
            for (var i = 0; i < x.length; i++) {
                sum += x[i];
            }
            ave = sum / x.length;
        }
        return ave;
    }

    /**
     * <p>
     * 由X值获得Y值
     * </p>
     * @param x
     *            当前X轴输入值,即为预测的月份
     * @param xx
     *            当前X轴输入值的前X数据点
     * @param a
     *            存储多项式系数的数组
     * @param m
     *            存储多项式的最高次数的数组
     *            <p>
     *            Date:2007-12-25 PM 20:07
     *            </p>
     * @return 对应X轴节点值的Y轴值
     */
    function getY(x, xx, a, m) {
        var y = 0;
        var ave = average(xx);

        var l = 0;
        for (var i = 0; i < m; i++) {
            l = a[0];
            if (i > 0) {
                y += a[i] * Math.pow((x - ave), i);
            }
        }
        return (y + l);
    }

    /**
     * 返回拟合后的点坐标数组
     * @param  {Array} arr 点坐标数组
     * @return {Array}     拟合后的点坐标数组
     */
    this.get = function(arr) {
        var arrX = [], arrY = [];

        for (var i = 0; i < arr.length; i++) {
            arrX.push(arr[i].x);
            arrY.push(arr[i].y);
        }

        var len = arrY.length;
        var m = 3;
        var a = new Array(arrX.length);
        var aa = PolyFit(arrX, arrY, len, a, m);
        var arrRes = [];
        for(var i = 0; i < len; i++){
           arrRes.push({x: arrX[i], y: getY(arrX[i], arrX, aa, m)});
        }

        return arrRes;
    };
}

// var arr = [{x: 10, y: 10},{x: 40, y: 90},{x: 70, y: 65},{x: 100, y: 15}];
// // 最小二乘法拟合
// var res = new PolyFitForcast().get(arr);
// var canvas = document.createElement('canvas');
// var ctx2d = canvas.getContext('2d');

// ctx2d.lineWidth = 1;

// ctx2d.strokeStyle = '#000';
// ctx2d.beginPath();
// ctx2d.moveTo(arr[0].x, arr[0].y);
// ctx2d.lineTo(arr[1].x, arr[1].y);
// ctx2d.lineTo(arr[2].x, arr[2].y);
// ctx2d.lineTo(arr[3].x, arr[3].y);
// ctx2d.stroke();

// ctx2d.strokeStyle = '#f00';
// ctx2d.beginPath();
// ctx2d.moveTo(res[0].x, res[0].y);
// ctx2d.lineTo(res[1].x, res[1].y);
// ctx2d.lineTo(res[2].x, res[2].y);
// ctx2d.lineTo(res[3].x, res[3].y);
// ctx2d.stroke();
// document.body.appendChild(canvas);

三.效果

把上面代码复制到浏览器console去掉注释再回车就能看到效果,当然,4个点的效果很差,下面给出10个点拟合的效果图:

curve

其实10个点拟合之后曲线还起来还不够圆滑,但如果有一个场景恰好能够提供很多个点,而且要求绘制速度还对精度要求不高的话,这种方式还是很不错的

四.更平滑自然的曲线拟合

最小二乘法是没指望了,除非反复进行多次拟合,但效果提升并不明显,曲线会越来越平,最后近似于直线,这不是我们想要的结果

还有一种方法叫贝塞尔曲线,除了一组信息点,还需要提供两个控制点,如果控制点选取得当就能实现相当完美的平滑曲线效果,难点是如何根据现有的信息点求出2个控制点,具体方法请查看Smoother Signatures

js实现曲线拟合(最小二乘法)》上有10条评论

    1. ayqy 文章作者

      根据返回的数组arrRes可以求,拟合出的曲线必定经过这些点,比如最后的例子是2次曲线,把几个点带入曲线方程一般式就可以解出来(2次曲线只需要3个点)

      回复
        1. ayqy 文章作者

          最小二乘法是根据已知的散点,求出一条把这些点一分为二的线,不要求穿过每一个点,如图

          线性回归的一种方法

          最小二乘就是根据误差最小求的(误差的平方和最小,可以参考基本原理

          拟合得到的点的数量和输入的点数量相同,因为y是根据输入的x得到的

          另外,刚才仔细看了下,发现曲线方程其实已经有了,就是getY函数(一个多项式),不用再把点带进去求了

          回复
      1. 辉仔

        我想请教下大神,你说的二次曲线的解析式是什么呀?可以发出来吗?代码没注释看着晕

        回复
        1. ayqy 文章作者

          方程:Y = a(0) + a(1) * (X – X1)+ a(2) * (X – X1)^2 + ….. .+ a(m) * (X – X1)^m

          注释中有的,this.get()中写死了默认var m = 3;,也就是三次曲线。其实不科学,应该根据待拟合点的特征,像几次m就是几。上面的m次多项式中,a是PolyFit()得到的系数数组,X1是x的平均值,m是最高次数

          回复
  1. 在线工具

    另外,在通过计算,给出最小二乘法的最小误差~~~就最好了~

    arrRes返回的时候拟合点,拟合点和被拟合的点数量是一样的吗?

    回复

发表评论

电子邮件地址不会被公开。 必填项已用*标注

*

code