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;
    }
    }
  • 相关阅读:
    SPOJ 694 (后缀数组) Distinct Substrings
    POJ 2774 (后缀数组 最长公共字串) Long Long Message
    POJ 3693 (后缀数组) Maximum repetition substring
    POJ 3261 (后缀数组 二分) Milk Patterns
    UVa 1149 (贪心) Bin Packing
    UVa 12206 (字符串哈希) Stammering Aliens
    UVa 11210 (DFS) Chinese Mahjong
    UVa (BFS) The Monocycle
    UVa 11624 (BFS) Fire!
    HDU 3032 (Nim博弈变形) Nim or not Nim?
  • 原文地址:https://www.cnblogs.com/sggggr/p/14247414.html
Copyright © 2011-2022 走看看