一、maven依赖

        <dependency>
            <groupId>org.apache.commons</groupId>
            <artifactId>commons-math3</artifactId>
            <version>3.6.1</version>
        </dependency>

二、线性拟合

math库已经实现了线性拟合，主要有3步

1、new 一个SimpleRegression 对象

2、向SimpleRegression 对象添加实验值addData

3、执行回归 regress，获取拟合值k,b和相关性系数R^2

说明一下R^2是衡量拟合度的一个参数，其值在0-1之间，值越接近1说拟合程度越好，

在医疗检测ESLA检测法中，实验要求R^2>0.98

        double[][] xy =new double[3][2];
        //TODO xy赋值
        SimpleRegression regression = new SimpleRegression();
        regression.addData(xy);
        //拟合参数
        RegressionResults results = regression.regress();
        double b = results.getParameterEstimate(0);
        double k = results.getParameterEstimate(1);
        //保留精度
        b = Double.parseDouble(DF.format(b));
        k = Double.parseDouble(DF.format(k));
        //相关性系数
        double rSquare = regression.getRSquare();

三、多项式拟合

       double[][] xy =new double[3][2];
        //TODO xy赋值
        //n 可变，两项式为2，三项式为3，以此类推，这只是设定的最大值，最终的多项式可能比这 
        //个值小
        PolynomialCurveFitter polynomialCurveFitter = PolynomialCurveFitter.create(n);
        WeightedObservedPoints weightedObservedPoints = new WeightedObservedPoints();
        for (int i = 0; i < xy.length; i++) {
            weightedObservedPoints.add(xy[i][0], xy[i][1]);
        }
        double[] fit = polynomialCurveFitter.fit(weightedObservedPoints.toList());
        StringBuilder polynomial = new StringBuilder("y = ");
        for (int i = fit.length - 1; i >= 0; i--) {
            fit[i] = Double.parseDouble(DF.format(fit[i]));
            if (fit[i] == 0) {
                continue;
            }
            // 多项式系数 数字
            String itemNum = fit[i] > 0 ? "+" : "";
            // 多项式 中的 x 次幂
            String itemXPower = i > 0 ? "*x^" + i : "";
            if (i == fit.length - 1) {
                // 最后一项常数
                polynomial.append(fit[i]).append(itemXPower);
            } else {
                // 带 x 的系数的项
                polynomial.append(itemNum).append(fit[i]).append(itemXPower);
            }
        }

四、四参数曲线拟合（自定义）

除了jdk自带的曲线拟合，jdk还提供了自定义的曲线拟合

1、实现ParametricUnivariateFunction 接口

value方法用于申明函数

gradient 则是对函数遍历参数（本例中是a、b、c、d）的求导

package cn.hwason.lab.common.math.function;

import org.apache.commons.math3.analysis.ParametricUnivariateFunction;

/**
 * f(x) = d + ((a - d) / (1 + (x / c) ^ b))
 * <p>
 * 指数函数的求导： f(x) = a^x => f'(x) = a^x * ln(a)
 * 密函数的求导：f(x) = x^a => f'(x) = a * x^(a-1)
 * 复合函数求导：u =f(x) v = g(u) => v' = g'(u) * u'
 * 分数的求导公式：f(x) = u(x) / v(x) => f'(x) = (u' * v - u * v') / v^2
 *
 * @author 印隽
 * @date 2023/7/24 9:43
 */

public class FourParameterFunction implements ParametricUnivariateFunction {

    /**
     * 定义公式
     *
     * @param x
     * @param parameters
     * @return
     */
    @Override
    public double value(double x, double... parameters) {
        double a = parameters[0];
        double b = parameters[1];
        double c = parameters[2];
        double d = parameters[3];
        return d + ((a - d) / (1 + Math.pow(x / c, b)));
    }

    /**
     * 对变量求导，需要和value()的定义顺序一致
     *
     * @param x
     * @param parameters
     * @return
     */
    @Override
    public double[] gradient(double x, double... parameters) {
        double a = parameters[0];
        double b = parameters[1];
        double c = parameters[2];
        double d = parameters[3];

        double[] gradients = new double[4];
        double den = 1 + Math.pow(x / c, b);

        gradients[0] = 1 / den; // 对 a 求导

        gradients[1] = -((a - d) * Math.pow(x / c, b) * Math.log(x / c)) / (den * den); // 对 b 求导

        gradients[2] = (b * Math.pow(x / c, b - 1) * (x / (c * c)) * (a - d)) / (den * den); // 对 c 求导

        gradients[3] = 1 - (1 / den); // 对 d 求导

        return gradients;

    }
}

2、接口调用

        double[][] xy = new double[5][2];
        //xy 赋值
        FourParameterFunction fourParameterFunction = new FourParameterFunction();
        double[] fit = CustomCurveFit.curveFit(xy, fourParameterFunction);
        //abcd和自己在value中定义的顺序是一致的
        double a = fit[0];
        double b = fit[1];
        double c = fit[2];
        double d = fit[3];
        //相关性系数
        double r_sqr = getR_sqr(xy, fourParameterFunction, a, b, c, d);

CustomCurveFit 定义

public class CustomCurveFit {

    public static double[] curveFit(double[][] xy, ParametricUnivariateFunction function) {
        double[] guess = new double[]{1, 1, 1, 1};
        SimpleCurveFitter curveFitter = SimpleCurveFitter.create(function, guess);
        WeightedObservedPoints observedPoints = new WeightedObservedPoints();
        for (int i = 0; i < xy.length; i++) {
            observedPoints.add(xy[i][0], xy[i][1]);
        }
        return curveFitter.fit(observedPoints.toList());
    }

}

3、相关性系数

自定义的相关性系数的计算方法，主要分3步

1.求 实验值与y的平均值 差 的平方和tss

2.求实验值与拟合的y值 差的平方和rss

3.R^2 = 1-tss/rss

  /**
     * 相关性 R^2 计算
     *
     * @param xy 用于拟合的实验值
     * @param parametricUnivariateFunction  自定义的拟合函数
     * @param params 拟合之后得到的fit[]值
     * @return
     */
    private static double getR_sqr(double[][] xy, ParametricUnivariateFunction parametricUnivariateFunction, double... params) {
        double sum_y = 0;
        for (int i = 0; i < xy.length; i++) {
            sum_y += xy[i][1];
        }
        //y的平均值
        double avg_y = sum_y / xy.length;

        //检测值（xy[i][1]）与y的平均值（avg_y）的 差值平方和
        double tss = 0;
        for (int i = 0; i < xy.length; i++) {
            tss += Math.pow(xy[i][1] - avg_y, 2);
        }
        //检测值（xy[i][1]）与拟合的y值（parametricUnivariateFunction.value(xy[i][0], params)）的差值平方和
        double rss = 0;
        for (int i = 0; i < xy.length; i++) {
            rss += Math.pow(xy[i][1] - parametricUnivariateFunction.value(xy[i][0], params), 2);
        }
        //R^2 计算
        return 1 - rss / tss;
    }

4、四参数的反函数

原函数 y = d + ((a - d) / (1 + (x / c) ^ b))

1. d移到等式左边  =>  y-d = (a-d)/(1+(x/c)^b)

2. a-d 移到等式左边 => (a-d)/(y-d) = 1+(x/c)^b

3. 1移到等式左边 => (a-d)/(y-d)-1 = (x/c)^b

4.等式两边同时取其对数函数 => In（ (a-d)/(y-d)-1）= In((x/c)^b)

5.利用对数函数的特性 loga^b = b*loga =>  In（ (a-d)/(y-d)-1）= b*In(x/c)

6. b移到等式左边 =>(1/b )*In（ (a-d)/(y-d)-1）= In(x/c)

7.两边同时取e的指数 => e^((1/b )*In（ (a-d)/(y-d)-1）) = e^( In(x/c))

8.利用公式e^Ina = a => e^((1/b )*In（ (a-d)/(y-d)-1）) = x/c

9. c移到等式左边 c*e^((1/b )*In（ (a-d)/(y-d)-1）)= x

10.这里已经得到反函数了x = c*e^((1/b )*In（ (a-d)/(y-d)-1）)，但是右边太复杂了，简化一下

11.利用公式e^Ina = a =>x= e^Inc * e^((1/b )*In（ (a-d)/(y-d)-1）)

12.利用公式(e^a)*(e^b) = e^(a+b) = > x = e^(lnc +(1/b )*In((a-d)/(y-d)-1))

13.分数加 => x = e^((b*lnc)/b +In((a-d)/(y-d)-1)/b)  =>  x = e^(((b*lnc) +In((a-d)/(y-d)-1))/b) 

14.利用公式 loga^b = b*loga => x = e^((Inc^b+In((a-d)/(y-d)-1))/b)

15.利用公式 loga+logb=logab => x=e^(In((c^b)*((a-d)/(y-d)-1)))/b)

16.简化一下 => x = e^ In(( (a-d)*c^b/(y-d) - c^b  )/b)

17.令  m = (a-d)*c^b ,n =  c^b => x= e^In((m/(y-d)-n)/b)

String.format("y=e^(In(%s/(x-%s)-%s)/%s)", Math.pow(c, b) * (a - d), d, Math.pow(c, b), b)

至此反函数就求解完成了，在ESLA检测法中很多试剂都是按照四参数进行曲线拟合，我们可以通过原函数 y = d + ((a - d) / (1 + (x / c) ^ b)) 拟合出a、b、c、d的值，得到反函数来展示拟合曲线