java实现离散点概率椭圆

package cn.headradio.supernova.monitoring.util;

import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.EigenDecomposition;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.stat.correlation.Covariance;

/**
 * @Author:Dava-J
 * @Description: 数组的一些算法
 * @Date: created in 2023年10月24日10:28:32
 * @Modified By:
 */
public class ArrayUtil {



    public static void main(String[] args) {
        // 输入数据
        double[][] data = {{1, 2}, {6, 2}, {8, 2}, {5, 7}, {1, 2}, {1, 2}, {1, 2}, {1, 2}, {7, 8}, {7, 7}};

        //Calculate the eigenvectors and eigenvalues
        // 转换为RealMatrix对象
        RealMatrix matrix = MatrixUtils.createRealMatrix(data);

        // 计算协方差矩阵
        Covariance covariance = new Covariance(matrix);
        RealMatrix covarianceMatrix = covariance.getCovarianceMatrix();

        // 打印结果
        for (int i = 0; i < covarianceMatrix.getRowDimension(); i++) {
            for (int j = 0; j < covarianceMatrix.getColumnDimension(); j++) {
                System.out.print(covarianceMatrix.getEntry(i, j) + " ");
            }
            System.out.println();
        }

        // 计算特征值和特征向量
        EigenDecomposition eigenDecomposition = new EigenDecomposition(covarianceMatrix);
        RealMatrix eigenvectors = verticalMatrixElements(eigenDecomposition.getV());
        RealMatrix eigenvalues = swapMatrixElements(eigenDecomposition.getD());

        // 打印结果
        System.out.println("特征值：");
        for (int i = 0; i < eigenvalues.getRowDimension(); i++) {

            System.out.println(eigenvalues.getEntry(i, i));
        }


        System.out.println("特征向量：");
        for (int i = 0; i < eigenvectors.getColumnDimension(); i++) {
            for (int j = 0; j < eigenvectors.getRowDimension(); j++) {
                System.out.print(eigenvectors.getEntry(j, i) + " ");
            }
            System.out.println();
        }


        double maxEigenvalue = max(eigenvalues);
        // 输出最大特征值及其索引
        System.out.println("最大特征值: " + maxEigenvalue);

        //Get the index of the largest eigenvector
        int[] indexs = new int[2];

        //获取最大值坐标
        for (int i = 0; i < eigenvalues.getRowDimension(); i++) {
            double[] row = eigenvalues.getRow(i);
            for (int j = 0; j < row.length; j++) {
                if (maxEigenvalue == eigenvalues.getEntry(i, j)) {
                    indexs[0] = i + 1;
                    indexs[1] = j + 1;
                }
            }

        }

        double[] largestEigenvec = eigenvectors.getColumn(indexs[0] - 1);
        //Get the smallest eigenvector and eigenvalue
        double smallestEigenval;
        double[] smallestEigenvec;
        if (indexs[0] == 1) {
            smallestEigenval = maxDouble(eigenvalues.getColumn(1));
            smallestEigenvec = eigenvectors.getRow(1);
        } else {
            smallestEigenval = maxDouble(eigenvalues.getColumn(0));
            smallestEigenvec = eigenvectors.getRow(0);
        }

        //This angle is between -pi and pi.
        //Let's shift it such that the angle is between 0 and 2pi
        double angle = Math.atan2(largestEigenvec[1], largestEigenvec[0]);

        if (angle < 0) {
            angle = angle + 2 * Math.PI;
        }

        //Get the coordinates of the data mean
        double xCount = 0;
        double yCount = 0;
        for (int i = 0; i < data.length; i++) {
            xCount += data[i][0];
            yCount += data[i][1];
        }
        double xAvg = xCount / data.length;
        double yAvg = yCount / data.length;

        double[] thetaGrid = getLinspace(0d, 2 * Math.PI);

        double X0 = 3.8;
        double Y0 = 3.6;

        double chisquareVal = 2.4477;
        double a = chisquareVal * Math.sqrt(maxEigenvalue);
        double b = chisquareVal * Math.sqrt(smallestEigenval);

        //the ellipse in x and y coordinates
        double[] ellipseXR = new double[thetaGrid.length];
        double[] ellipseYR = new double[thetaGrid.length];
        for (int i = 0; i < thetaGrid.length; i++) {
            ellipseXR[i] = a * Math.cos(thetaGrid[i]);
        }
        for (int i = 0; i < thetaGrid.length; i++) {
            ellipseYR[i] = b * Math.sin(thetaGrid[i]);
        }

        //Define a rotation matrix
        double[][] r = new double[2][2];
        r[0][0] = Math.cos(angle);
        r[0][1] = Math.sin(angle);
        r[1][0] = -Math.sin(angle);
        r[1][1] = Math.cos(angle);

        double[][] merge = new double[thetaGrid.length][2];
        for (int i = 0; i < merge.length; i++) {
            merge[i][0] = ellipseXR[i];
            merge[i][1] = ellipseYR[i];
        }

        double[][] result = multiplyMatrices(merge, r);
        //let's rotate the ellipse to some angle phi


        System.out.println();

    }


    /**
     * 矩阵相乘
     * @param matrix1
     * @param matrix2
     * @return
     */
    public static double[][] multiplyMatrices(double[][] matrix1, double[][] matrix2) {
        int rows1 = matrix1.length;
        int cols1 = matrix1[0].length;
        int cols2 = matrix2[0].length;

        double[][] result = new double[rows1][cols2];

        for (int i = 0; i < rows1; i++) {
            for (int j = 0; j < cols2; j++) {
                for (int k = 0; k < cols1; k++) {
                    result[i][j] += matrix1[i][k] * matrix2[k][j];
                }
            }
        }

        return result;
    }


    /**
     * 获取
     * @return
     */
    public static double[] getLinspace(double startNum, double intervalSize){
        int size = 100;
        intervalSize = intervalSize * 0.01;
        double[] result = new double[size];
        for (int i = 0; i < size; i++) {
            result[i] = startNum + i * intervalSize;
        }
        double[] lin = new double[]{0,0.0634665182543393,0.126933036508679,0.190399554763018,0.253866073017357,0.317332591271696,0.380799109526036,0.444265627780375,0.507732146034714,0.571198664289053,0.634665182543393,0.698131700797732,0.761598219052071,0.825064737306410,0.888531255560750,0.951997773815089,
                1.01546429206943,1.07893081032377,1.14239732857811,1.20586384683245,1.26933036508679,1.33279688334112,1.39626340159546,1.45972991984980,1.52319643810414,1.58666295635848,1.65012947461282,1.71359599286716,1.77706251112150,1.84052902937584,1.90399554763018,1.96746206588452,
                2.03092858413886,2.09439510239320,2.15786162064753,2.22132813890187,2.28479465715621,2.34826117541055,2.41172769366489,2.47519421191923,2.53866073017357,2.60212724842791,2.66559376668225,2.72906028493659,2.79252680319093,2.85599332144527,2.91945983969961,2.98292635795395,
                3.04639287620828,3.10985939446262,3.17332591271696,3.23679243097130,3.30025894922564,3.36372546747998,3.42719198573432,3.49065850398866,3.55412502224300,3.61759154049734,3.68105805875168,3.74452457700602,3.80799109526036,3.87145761351469,3.93492413176903,3.99839065002337,
                4.06185716827771,4.12532368653205,4.18879020478639,4.25225672304073,4.31572324129507,4.37918975954941,4.44265627780375,4.50612279605809,4.56958931431243,4.63305583256677,4.69652235082111,4.75998886907544,4.82345538732978,4.88692190558412,4.95038842383846,
                5.01385494209280,5.07732146034714,5.14078797860148,5.20425449685582,5.26772101511016,5.33118753336450,5.39465405161884,5.45812056987318,5.52158708812751,5.58505360638185,5.64852012463619,5.71198664289053,5.77545316114487,5.83891967939921,5.90238619765355,5.96585271590789,
                6.02931923416223,6.09278575241657,6.15625227067091,6.21971878892525,6.28318530717959};
        return lin;
    }
    /**
     * 获取数组中最大的值
     *
     * @param arr
     * @return
     */
    public static double maxDouble(double[] arr) {
        double max = arr[0];
        for (int i = 0; i < arr.length; i++) {
            if (arr[i] > max) {
                max = arr[i];
            }
        }
        return max;
    }


    /**
     * 取最大值
     *
     * @param eigenvalues
     * @return
     */
    public static double max(RealMatrix eigenvalues) {
        // 使用EigenDecomposition对矩阵进行特征值分解
        EigenDecomposition eigenDecomposition1 = new EigenDecomposition(eigenvalues);

        // 获取特征值
        double[] eigenvalues1 = eigenDecomposition1.getRealEigenvalues();

        // 找到最大特征值和其索引
        double maxEigenvalue = Double.NEGATIVE_INFINITY;

        for (int i = 0; i < eigenvalues1.length; i++) {
            if (eigenvalues1[i] > maxEigenvalue) {
                maxEigenvalue = eigenvalues1[i];
            }
        }
        return maxEigenvalue;
    }

    /**
     * @param originalMatrix
     * @return
     */
    public static RealMatrix swapMatrixElements(RealMatrix originalMatrix) {
        double x1 = originalMatrix.getEntry(0, 0);
        double y1 = originalMatrix.getEntry(0, 1);
        double x2 = originalMatrix.getEntry(1, 0);
        double y2 = originalMatrix.getEntry(1, 1);

        RealMatrix swappedMatrix = new Array2DRowRealMatrix(new double[][]{
                {y2, x2},
                {y1, x1}
        });

        return swappedMatrix;
    }

    /**
     * @param originalMatrix
     * @return
     */
    public static RealMatrix verticalMatrixElements(RealMatrix originalMatrix) {
        double x1 = originalMatrix.getEntry(0, 0);
        double y1 = originalMatrix.getEntry(0, 1);
        double x2 = originalMatrix.getEntry(1, 0);
        double y2 = originalMatrix.getEntry(1, 1);

        RealMatrix swappedMatrix = new Array2DRowRealMatrix(new double[][]{
                {x2, y2},
                {-x1, -y1}
        });

        return swappedMatrix;
    }

    // 打印矩阵
    public static void printMatrix(RealMatrix matrix) {
        for (double[] row : matrix.getData()) {
            for (double value : row) {
                System.out.print(value + " ");
            }
            System.out.println();
        }
        System.out.println();
    }
}