zoukankan      html  css  js  c++  java
  • 各种拟合,一元、多元、对数、指数、单峰、自定义拟合【转】

    其中引用到了apache的common-math的jar包,主要用于矩阵运算,下载地址:
    
    http://commons.apache.org/proper/commons-math/userguide/fitting.html
    
    import org.apache.commons.math3.fitting.PolynomialCurveFitter;
    import org.apache.commons.math3.fitting.WeightedObservedPoints;
    import org.apache.commons.math3.linear.ArrayRealVector;
    import org.apache.commons.math3.linear.DecompositionSolver;
    import org.apache.commons.math3.linear.LUDecomposition;
    import org.apache.commons.math3.linear.MatrixUtils;
    import org.apache.commons.math3.linear.RealMatrix;
    import org.apache.commons.math3.linear.RealVector;
    
    
    public class MathUtil 
    {
    /**
    * 一元线性拟合
    * y = a + b*x
    * @param x  
    * @param y
    * reuslt[0] = a
    * result[1] = b
    * result[2] 相关系数
    * result[3] 决定系数
    * result[4] 点数量(长度)
    * result[5] 自由度
    */
    public static double[] lineFitting(double x[], double y[])
    {
    int size = x.length;
    double xmean = 0.0;
    double ymean = 0.0;
    double xNum = 0.0;
    double yNum = 0.0;
    double xyNum = 0;
    double xNum2 = 0;
    double yNum2 = 0;
    double rss = 0;
    double tss = 0;
    double result[] = new double[6];
    
    for (int i = 0; i < size; i++)
    {
    xmean += x[i];
    ymean += y[i];
    xNum2 += x[i] * x[i];
    yNum2 += y[i] * y[i];
    xyNum += x[i] * y[i];
    }
    xNum = xmean;
    yNum = ymean;
    xmean /= size;
    ymean /= size;
    
    double sumx2 = 0.0f;
    double sumxy = 0.0f;
    for (int i = 0; i < size; i++)
    {
    sumx2 += (x[i] - xmean) * (x[i] - xmean);
    sumxy += (y[i] - ymean) * (x[i] - xmean);
    }
    
    
    double b = sumxy / sumx2;
    double a = ymean - b * xmean;
    
    
    result[0] = a;
    result[1] = b;
    System.out.println("a = " + a + ", b=" + b);
    
    double correlation = (xyNum - xNum * yNum / size) / Math.sqrt((xNum2 - xNum * xNum / size) * (yNum2 - yNum * yNum / size));
    
    System.out.println("相关系数:" + correlation);
    result[2] = correlation;
    
    for (int i = 0; i < size; i++)
    {
    rss += (y[i] - (a + b * x[i])) * (y[i] - (a + b * x[i]));
    tss += (y[i] - ymean) * (y[i] - ymean);
    }
    
    double r2 = 1 - (rss / (size - 1 - 1)) / (tss/ (size - 1));
    
    result[3] = r2;
    System.out.println("决定系数" + r2);
    
    result[4] = x.length;
    result[5] = x.length - 1 - 1;
    
    return result;
    }
    
    /**
    * 多元线性拟合
    * y = a + b*x1 + c*x2
    * @param x  
    * @param y
    * result[0] = a
    * result[b] = b
    * .
    * .
    * .
    * result[len - 4] 点数
    * result[len - 3] 自由度
    * result[len - 2] 残差平方和
    * result[len - 1] 确定系数
    */
    public static double[] lineFitting2(double x[][], double y[])
    {
    double[] a = new double[x.length + 1];
    double[] v = new double[2];  // 这里的2为m
    double[] dt = new double[4];
    
    
    line2sqt(x, y, 2, 11, a, dt, v);
    
    int i;
    
    System.out.println("残差平方和:" + dt[0]);
    
    double temp = a[a.length - 1];
    // 更换输出位置,把常数放到第一位
    for(i = a.length - 1; i > 0; i--)
    {
    a[i] = a[i - 1];
    }
    
    a[0] = temp;
    
    double[] result = new double[x.length + 5];
    
    for (i = 0; i <= x.length; i++)
    {
    result[i] = a[i];
    }
    
    result[x.length + 1] = y.length;
    result[x.length + 2] = y.length - x.length;
    result[x.length + 3] = dt[0];
    result[x.length + 4] = getLine2R(x, y, a, x.length);
    
    System.out.println("校正确定系数:" + result[x.length + 4]);
    
    return result;
    }
    
    /**
    * 多项式拟合
    * y = a + b*x1 + c*x1^2……
    * result[0] = a
    * result[1] = b
    * .
    * .
    * .
    * result[n + 1] 点数
    * result[n + 2] 自由度
    * result[n + 3] 确定系数
    * @param n 几级
    * @return
    */
    public static double[] dxsFitting(double x[], double y[], int n)
    {
    double result[] = new double[n + 4];
    
    WeightedObservedPoints obs = new WeightedObservedPoints();
    
    for(int i = 0; i < x.length; i++)
    {
    obs.add(x[i], y[i]);
    }
    
    // Instantiate a third-degree polynomial fitterm.
    final PolynomialCurveFitter fitter = PolynomialCurveFitter.create(n);
    
    
    // Retrieve fitted parameters (coefficients of the polynomial function).
    final double[] coeff = fitter.fit(obs.toList());
    
    for(int i = 0; i < coeff.length; i++)
    {
    result[i] = coeff[i];
    }
    
    double s = getDxsR(x, y, coeff, 5);
    
    System.out.println("确定系数:" + s);
    
    result[n + 1] = x.length;
    result[n + 2] = x.length - n - 1;
    result[n + 3] = s;
    
    return result;
    }
    
    /**
    * 指数拟合
    * y = b*exp(ax)
    * result[0] = a
    * result[1] = b
    * result[2] 数据点数
    * result[3] 自由度
    * result[4] 确定系数
    * @param x
    * @param y
    * @return
    */
    public static double[] expFitting(double x[], double y[])
    {
    int size = x.length;
    double xmean = 0.0;
    double ymean = 0.0;
    double rss = 0;
    double tss = 0;
    double result[] = new double[5];
    
    for (int i = 0; i < size; i++)
    {
    xmean += x[i];
    y[i] = Math.log(y[i]);
    ymean += y[i];
    }
    xmean /= size;
    ymean /= size;
    
    double sumx2 = 0.0f;
    double sumxy = 0.0f;
    for (int i = 0; i < size; i++)
    {
    sumx2 += (x[i] - xmean) * (x[i] - xmean);
    sumxy += (y[i] - ymean) * (x[i] - xmean);
    }
    
    
    double b = sumxy / sumx2;
    double a = ymean - b * xmean;
    
    for (int i = 0; i < size; i++)
    {
    rss += (y[i] - (a + b * x[i])) * (y[i] - (a + b * x[i]));
    tss += (y[i] - ymean) * (y[i] - ymean);
    }
    
    double r2 = 1 - (rss / (size - 1 - 1)) / (tss/ (size - 1));
    
    System.out.println("决定系数" + r2);
    
    
    a = Math.exp(a);
    
    System.out.println("a = " + a + ";b= " + b);
    
    result[0] = a;
    result[1] = b;
    result[2] = x.length;
    result[3] = x.length - 2;
    result[4] = r2;
    
    return result;
    }
    
    /**
    * 对数拟合
    * y = ln(a * x + b)
    * result[0] = a
    * result[1] = b
    * result[2] 数据点数
    * result[3] 自由度
    * result[4] 确定系数
    * @param x
    * @param y
    */
    public static double[] logFitting(double x[], double y[])
    {
    int size = x.length;
    double xmean = 0.0;
    double ymean = 0.0;
    double rss = 0;
    double tss = 0;
    double result[] = new double[5];
    
    for (int i = 0; i < size; i++)
    {
    xmean += x[i];
    y[i] = Math.exp(y[i]);
    ymean += y[i];
    }
    xmean /= size;
    ymean /= size;
    
    double sumx2 = 0.0f;
    double sumxy = 0.0f;
    for (int i = 0; i < size; i++)
    {
    sumx2 += (x[i] - xmean) * (x[i] - xmean);
    sumxy += (y[i] - ymean) * (x[i] - xmean);
    }
    
    
    double b = sumxy / sumx2;
    double a = ymean - b * xmean;
    
    for (int i = 0; i < size; i++)
    {
    rss += (y[i] - (b + a * x[i])) * (y[i] - (b + a * x[i]));
    tss += (y[i] - ymean) * (y[i] - ymean);
    }
    
    double r2 = 1 - (rss / (size - 1 - 1)) / (tss/ (size - 1));
    
    System.out.println("决定系数" + r2);
    
    
    System.out.println("a = " + a + ";b= " + b);
    
    result[0] = a;
    result[1] = b;
    result[2] = x.length;
    result[3] = x.length - 2;
    result[4] = r2;
    
    return result;
    
    }
    
    /**
    * 参考网址:http://wenku.baidu.com/link?url=T-xy9CUR_hVlCpA5o6xjqWb-Z1o5jraGEg-OHEtDgFLEynf6HedV68wMwC5u7RTs7PN5kiYq3qjyiddMveDFywLiuXai11MisNrNvv4KjAW&qq-pf-to=pcqq.c2c
    * 峰拟合
    * y = y(max) * exp[-(x - x(max))^2/S]
    * @param x
    * @param y
    * result[0] = x(max)
    * result[1] = y(max)
    * result[2] = S
    * result[3] 点数
    * result[4] 自由度
    * result[5] 确定系数
    */
    public static double[] peakFitting(double x[], double y[])
    {
    int size = x.length;
    double maxX = 0;
    double maxY = 0;
    double minY = Integer.MAX_VALUE;
    double result[] = new double[6];
    double[][] left = new double[x.length][3];
    double[][] right = new double[x.length][1];
    
    for (int i = 0; i < size; i++)
    {
                if(y[i] > maxY)
                {
                maxX = x[i];
                    maxY = y[i];
                }
    
    
                if(y[i] < minY)
                {
                    minY = y[i];
                }
                
                for(int j = 0; j < 3; j++)
    {
    left[i][j] = getPeakValue(j, x[i]);
    }
    
    right[i][0] = Math.log(y[i]);
    }
    
    RealMatrix leftMatrix = MatrixUtils.createRealMatrix(left);
    RealMatrix rightMatrix = MatrixUtils.createRealMatrix(right);
    RealMatrix leftMatrix1 = leftMatrix.transpose();
    RealMatrix m = leftMatrix1.multiply(leftMatrix);
    
    RealMatrix tMatrix = new LUDecomposition(m).getSolver().getInverse().multiply(leftMatrix1).multiply(rightMatrix);
    
    result[0] = maxX;
    result[1] = maxY;
    result[2] = - 1.0 / tMatrix.getEntry(2, 0);
    
    System.out.println(result[0] + " " + result[1] + " " + result[2]);
    double aaaa = getPearR(x, y, result, 2);
    System.out.println("确定系数:" + aaaa);
    
    result[3] = x.length;
    result[4] = x.length - 1 - 2;
    result[5] = aaaa;
    
    return result;
    }
    
    /**
    * 用户自定义拟合
    * @param x
    * @param y
    * @param sf   算法
    *  <input type="checkbox" value="0" />常数
    <input type="checkbox" value="1" />x
    <input type="checkbox" value="2" />x^2
    <input type="checkbox" value="3" />x^3
    <input type="checkbox" value="4" />exp(x)
    <input type="checkbox" value="5" />ln(x)
    <input type="checkbox" value="6" />sin(x)
    <input type="checkbox" value="7" />cos(x)
    * 
    * result[0] = 第一项系数
    * result[1] = 第二项系数
    * .
    * .
    * .
    * result[n] = 第N项系数
    * result[sf.length]   点数
    * result[sf.length+1] 自由度
    * result[sf.length+2] 确定系数
    * @return
    */
    public static double[] userDefineFitting(double x[], double y[], int sf[])
    {
    double[][] left = new double[sf.length][sf.length];
    double[] right = new double[sf.length];
    double[] result = new double[sf.length + 3];
    
    result[sf.length] = x.length;
    boolean containZero = false;
    
    for(int i = 0; i < sf.length; i++)
    {
    double yValue = 0;
    
    // 数值中包括0
    if(sf[i] == 0)
    {
    containZero = true;
    }
    
    for(int j = 0; j < sf.length; j++)
    {
    double xValue = 0;
    
    // 计算X的的司格马
    for(int k = 0; k < x.length; k++)
    {
    xValue += getUserDefineValue(sf[i], x[k]) * getUserDefineValue(sf[j], x[k]);
    }
    
    left[i][j] = xValue;
    }
    
    // 计算Y的的司格马
    for(int k = 0; k < x.length; k++)
    {
    yValue += y[k] * getUserDefineValue(sf[i], x[k]);
    }
    
    right[i] = yValue;
    }
    
    // 计算自由度
    result[sf.length + 1] = x.length - sf.length - 1;
    
    // 计算自由度
    if(containZero)
    {
    result[sf.length + 1] = x.length - sf.length;
    }
    
    // 矩阵求解
    RealMatrix leftMatrix = MatrixUtils.createRealMatrix(left);
    DecompositionSolver solver = new LUDecomposition(leftMatrix).getSolver();
    RealVector constants = new ArrayRealVector(right, false);
    
    RealVector solution = solver.solve(constants);
    
    System.out.println(solution);
    
    // 获得系数值
    for(int i = 0; i < solution.getDimension(); i++)
    {
    result[i] = solution.getEntry(i);
    }
    
    double rss = 0, tss = 0, ymean = 0;
    
    for (int i = 0; i < x.length; i++)
    {
    ymean += y[i];
    }
    
    ymean /= x.length;
    
    for (int i = 0; i < x.length; i++)
    {
    rss += (y[i] - getUserDefineValueByX(sf, x[i], solution)) * (y[i] - getUserDefineValueByX(sf, x[i], solution));
    tss += (y[i] - ymean) * (y[i] - ymean);
    }
    
    double r2 = 1 - (rss / result[sf.length + 1]) / (tss / (x.length - 1));
    
    System.out.println("决定系数" + r2);
    
    result[sf.length + 2] = r2;
    
    return result;
    }
    
    /**
    * 用户自定义函数,传入X值获得Y值
    * @param n
    * @param x
    * @param solution
    * @return
    */
    private static double getUserDefineValueByX(int[] n, double x, RealVector solution)
    {
    double value = 0;
    
    for(int i = 0; i < n.length; i++)
    {
    value += getUserDefineValue(n[i], x) * solution.getEntry(i);
    }
    
    return value;
    }
    
    /**
    * 获得用户自定义的值, 
    * <input type="checkbox" value="0" />常数
    <input type="checkbox" value="1" />x
    <input type="checkbox" value="2" />x^2
    <input type="checkbox" value="3" />x^3
    <input type="checkbox" value="4" />exp(x)
    <input type="checkbox" value="5" />ln(x)
    <input type="checkbox" value="6" />sin(x)
    <input type="checkbox" value="7" />cos(x)
    * @param i
    * @param x
    * @return
    */
    private static double getUserDefineValue(int i, double x)
    {
    // 常数
    if(i == 0)
    {
    return 1;
    }else if(i == 1)
    {
    // x
    return x;
    }else if(i == 2)
    {
    // x^2
    return x * x;
    }else if(i == 3)
    {
    // x^3
    return x * x * x;
    }else if(i == 4)
    {
    //exp(x)
    return Math.pow(Math.E, x);
    }else if(i == 5)
    {
    // ln(x)
    return Math.log(x);
    }else if(i == 6)
    {
    // sin(x)
    return Math.sin(x);
    }else if(i == 7)
    {
    // cos(x)
    return Math.cos(x);
    }
    
    return 0;
    }
    
    /**
    * 获得决定系数
    * @param coeff
    * @return
    */
    private static double getPearR(double x[], double y[], double coeff[], int n)
    {
    int size = x.length;
    double ymean = 0.0;
    double rss = 0;
    double tss = 0;
    
    for (int i = 0; i < size; i++)
    {
    ymean += y[i];
    }
    
    ymean /= size;
    
    for (int i = 0; i < size; i++)
    {
    rss += (y[i] - getPeakValueByX(x[i], coeff)) * (y[i] - getPeakValueByX(x[i], coeff));
    tss += (y[i] - ymean) * (y[i] - ymean);
    }
    
    double r2 = 1 - (rss / (size - n - 1)) / (tss / (size - 1));
    
    return r2;
    }
    
    /**
    * 通过X获得Y的值
    * y = y(max) * exp[-(x - x(max))^2/S]
    * @param x
    * @param y
    * result[0] = x(max)
    * result[1] = y(max)
    * result[2] = S
    * @return
    */
    private static double getPeakValueByX(double x, double result[])
    {
    double b = Math.exp(-Math.pow((x - result[0]), 2) / result[2]) * result[1];
    
    return b;
    }
    
    private static double getPeakValue(int i, double x)
    {
    // 常数
    if(i == 0)
    {
    return 1;
    }else if(i == 1)
    {
    return x;
    }else if(i == 2)
    {
    // x^2
    return x * x;
    }
    
    return 0;
    }
    
    /**
    * 获得决定系数
    * @param coeff
    * @return
    */
    private static double getDxsR(double x[], double y[], double coeff[], int n)
    {
    int size = x.length;
    double ymean = 0.0;
    double rss = 0;
    double tss = 0;
    
    for (int i = 0; i < size; i++)
    {
    ymean += y[i];
    }
    
    ymean /= size;
    
    for (int i = 0; i < size; i++)
    {
    rss += (y[i] - getDxsValueByX(x[i], coeff)) * (y[i] - getDxsValueByX(x[i], coeff));
    tss += (y[i] - ymean) * (y[i] - ymean);
    }
    
    double r2 = 1 - (rss / (size - n - 1)) / (tss / (size - 1));
    
    return r2;
    }
    
    /**
    * 通过X获得Y的值
    * @param x
    * @param coeff
    * @return
    */
    private static double getDxsValueByX(double x, double coeff[])
    {
    int size = coeff.length;
    double result = coeff[0];
    
    for(int i = 1; i < size; i++)
    {
    result += coeff[i] * Math.pow(x, i);
    }
    
    return result;
    }
    
    
    /**
    * 多元线性回归分析
    * 
    * @param x[m][n]
    *            每一列存放m个自变量的观察值
    * @param y[n]
    *            存放随即变量y的n个观察值
    * @param m
    *            自变量的个数
    * @param n
    *            观察数据的组数
    * @param a
    *            返回回归系数a0,...,am
    * @param dt[4]
    *            dt[0]偏差平方和q,dt[1] 平均标准偏差s dt[2]返回复相关系数r dt[3]返回回归平方和u
    * @param v[m]
    *            返回m个自变量的偏相关系数
    */
    private static void line2sqt(double[][] x, double[] y, int m, int n, double[] a, double[] dt, double[] v)
    {
    int i, j, k, mm;
    double q, e, u, p, yy, s, r, pp;
    double[] b = new double[(m + 1) * (m + 1)];
    mm = m + 1;
    b[mm * mm - 1] = n;
    for (j = 0; j <= m - 1; j++)
    {
    p = 0.0;
    for (i = 0; i <= n - 1; i++)
    p = p + x[j][i];
    b[m * mm + j] = p;
    b[j * mm + m] = p;
    }
    for (i = 0; i <= m - 1; i++)
    for (j = i; j <= m - 1; j++)
    {
    p = 0.0;
    for (k = 0; k <= n - 1; k++)
    p = p + x[i][k] * x[j][k];
    b[j * mm + i] = p;
    b[i * mm + j] = p;
    }
    a[m] = 0.0;
    for (i = 0; i <= n - 1; i++)
    a[m] = a[m] + y[i];
    for (i = 0; i <= m - 1; i++)
    {
    a[i] = 0.0;
    for (j = 0; j <= n - 1; j++)
    a[i] = a[i] + x[i][j] * y[j];
    }
    
    chlk(b, mm, 1, a);
    
    yy = 0.0;
    for (i = 0; i <= n - 1; i++)
    yy = yy + y[i] / n;
    q = 0.0;
    e = 0.0;
    u = 0.0;
    for (i = 0; i <= n - 1; i++)
    {
    p = a[m];
    for (j = 0; j <= m - 1; j++)
    p = p + a[j] * x[j][i];
    q = q + (y[i] - p) * (y[i] - p);
    e = e + (y[i] - yy) * (y[i] - yy);
    u = u + (yy - p) * (yy - p);
    }
    s = Math.sqrt(q / n);
    r = Math.sqrt(1.0 - q / e);
    for (j = 0; j <= m - 1; j++)
    {
    p = 0.0;
    for (i = 0; i <= n - 1; i++)
    {
    pp = a[m];
    for (k = 0; k <= m - 1; k++)
    if (k != j)
    pp = pp + a[k] * x[k][i];
    p = p + (y[i] - pp) * (y[i] - pp);
    }
    v[j] = Math.sqrt(1.0 - q / p);
    }
    dt[0] = q;
    dt[1] = s;
    dt[2] = r;
    dt[3] = u;
    }
    
    
    private static int chlk(double[] a, int n, int m, double[] d)
    {
    int i, j, k, u, v;
    if ((a[0] + 1.0 == 1.0) || (a[0] < 0.0))
    {
    System.out.println("fail
    ");
    return (-2);
    }
    a[0] = Math.sqrt(a[0]);
    for (j = 1; j <= n - 1; j++)
    a[j] = a[j] / a[0];
    for (i = 1; i <= n - 1; i++)
    {
    u = i * n + i;
    for (j = 1; j <= i; j++)
    {
    v = (j - 1) * n + i;
    a[u] = a[u] - a[v] * a[v];
    }
    if ((a[u] + 1.0 == 1.0) || (a[u] < 0.0))
    {
    System.out.println("fail
    ");
    return (-2);
    }
    a[u] = Math.sqrt(a[u]);
    if (i != (n - 1))
    {
    for (j = i + 1; j <= n - 1; j++)
    {
    v = i * n + j;
    for (k = 1; k <= i; k++)
    a[v] = a[v] - a[(k - 1) * n + i] * a[(k - 1) * n + j];
    a[v] = a[v] / a[u];
    }
    }
    }
    for (j = 0; j <= m - 1; j++)
    {
    d[j] = d[j] / a[0];
    for (i = 1; i <= n - 1; i++)
    {
    u = i * n + i;
    v = i * m + j;
    for (k = 1; k <= i; k++)
    d[v] = d[v] - a[(k - 1) * n + i] * d[(k - 1) * m + j];
    d[v] = d[v] / a[u];
    }
    }
    for (j = 0; j <= m - 1; j++)
    {
    u = (n - 1) * m + j;
    d[u] = d[u] / a[n * n - 1];
    for (k = n - 1; k >= 1; k--)
    {
    u = (k - 1) * m + j;
    for (i = k; i <= n - 1; i++)
    {
    v = (k - 1) * n + i;
    d[u] = d[u] - a[v] * d[i * m + j];
    }
    v = (k - 1) * n + k - 1;
    d[u] = d[u] / a[v];
    }
    }
    return (2);
    }
    
    /**
    * 获得相关系数
    * @param coeff
    * @return
    */
    private static double getLine2R(double x[][], double y[], double a[], int n)
    {
    int size = x[0].length;
    double ymean = 0.0;
    double rss = 0;
    double tss = 0;
    
    for (int i = 0; i < size; i++)
    {
    ymean += y[i];
    }
    
    ymean /= size;
    
    for (int i = 0; i < size; i++)
    {
    rss += (y[i] - getLine2ValueByX(x, a, i)) * (y[i] - getLine2ValueByX(x, a, i));
    tss += (y[i] - ymean) * (y[i] - ymean);
    }
    
    double r2 = 1 - (rss / (size - n - 1)) / (tss / (size - 1));
    
    return r2;
    }
    
    /**
    * 通过X获得Y的值
    * @param x
    * @return
    */
    private static double getLine2ValueByX(double x[][], double a[], int n)
    {
    int size = a.length;
    double result = a[0];
    
    for(int i = 1; i < size; i++)
    {
    result += x[i - 1][n] * a[i];
    }
    
    return result;
    }
    }
  • 相关阅读:
    RabbitMQ教程
    设置文件夹共享及所有访问权限
    一键操作IIS
    C#操作IIS服务
    【01】浅谈HTTP在WebApi开发中的运用
    Socket通信原理
    C# 通过JObject解析json对象
    MVC小知识整理
    Global.asax文件的作用
    iOS开发日记9-终端命令
  • 原文地址:https://www.cnblogs.com/sggggr/p/14247414.html
Copyright © 2011-2022 走看看