java实现GMM算法
public class Parameter {private ArrayList<ArrayList<Double>> pMiu; // 均值参数k个分布的中心点，每个中心点d维private ArrayList<Double> pPi; // k个GMM的权值private ArrayList<ArrayList<ArrayList<Double>>> pSigma; // k类GMM的协方差矩阵,d*d*kpublic ArrayList<ArrayList<Double>> getpMiu() {return pMiu;}public void setpMiu(ArrayList<ArrayList<Double>> pMiu) {this.pMiu = pMiu;}public ArrayList<Double> getpPi() {return pPi;}public void setpPi(ArrayList<Double> pPi) {this.pPi = pPi;}public ArrayList<ArrayList<ArrayList<Double>>> getpSigma() {return pSigma;}public void setpSigma(ArrayList<ArrayList<ArrayList<Double>>> pSigma) {this.pSigma = pSigma;}}

?

核心代码如下：

public class GMMAlgorithm {/** * * @Title: GMMCluster * @Description: GMM聚类算法的实现类，返回每条数据的类别(0~k-1)* @return int[]* @throws */public int[] GMMCluster(ArrayList<ArrayList<Double>>dataSet, ArrayList<ArrayList<Double>> pMiu, int dataNum, int k, int dataDimen) {Parameter parameter = iniParameters(dataSet, dataNum, k, dataDimen);double Lpre = -1000000; // 上一次聚类的误差double threshold = 0.0001;while(true) {ArrayList<ArrayList<Double>> px = computeProbablity(dataSet, pMiu, dataNum, k, dataDimen);double[][] pGama = new double[dataNum][k];for(int i = 0; i < dataNum; i++) {for(int j = 0; j < k; j++) {pGama[i][j] = px.get(i).get(j) * parameter.getpPi().get(j);}}double[] sumPGama = GMMUtil.matrixSum(pGama, 2);for(int i = 0; i < dataNum; i++) {for(int j = 0; j < k; j++) {pGama[i][j] = pGama[i][j] / sumPGama[i];}}double[] NK = GMMUtil.matrixSum(pGama, 1); // 第k个高斯生成每个样本的概率的和，所有Nk的总和为N// 更新pMiudouble[] NKReciprocal = new double[NK.length];for(int i = 0; i < NK.length; i++) {NKReciprocal[i] = 1 / NK[i];}double[][] pMiuTmp = GMMUtil.matrixMultiply(GMMUtil.matrixMultiply(GMMUtil.diag(NKReciprocal), GMMUtil.matrixReverse(pGama)), GMMUtil.toArray(dataSet));// 更新pPiedouble[][] pPie = new double[k][1];for(int i = 0; i < NK.length; i++) {pPie[i][1] = NK[i] / dataNum;}// 更新k个pSigmadouble[][][] pSigmaTmp = new double[dataDimen][dataDimen][k];for(int i = 0; i < k; i++) {double[][] xShift = new double[dataNum][dataDimen];for(int j = 0; j < dataNum; j++) {for(int l = 0; l < dataDimen; l++) {xShift[j][l] = pMiuTmp[i][l];}}double[] pGamaK = new double[dataNum]; // 第k条pGama值for(int j = 0; j < dataNum; j++) {pGamaK[j] = pGama[j][i];}double[][] diagPGamaK = GMMUtil.diag(pGamaK);double[][] pSigmaK = GMMUtil.matrixMultiply(GMMUtil.matrixReverse(xShift), (GMMUtil.matrixMultiply(diagPGamaK, xShift)));for(int j = 0; j < dataDimen; j++) {for(int l = 0; l < dataDimen; l++) {pSigmaTmp[j][l][k] = pSigmaK[j][l] / NK[i];}}}// 判断是否迭代结束double[][] a = GMMUtil.matrixMultiply(GMMUtil.toArray(px), pPie);for(int i = 0; i < dataNum; i++) {a[i][0] = Math.log(a[i][0]);}double L = GMMUtil.matrixSum(a, 1)[0];if(L - Lpre < threshold) {break;}Lpre = L;}return null;}/** * * @Title: computeProbablity * @Description: 计算每个节点（共n个）属于每个分布（k个）的概率* @return ArrayList<ArrayList<Double>>* @throws */public ArrayList<ArrayList<Double>> computeProbablity(ArrayList<ArrayList<Double>>dataSet, ArrayList<ArrayList<Double>> pMiu, int dataNum, int k, int dataDimen) {double[][] px = new double[dataNum][k]; // 每条数据属于每个分布的概率 int[] type = getTypes(dataSet, pMiu, k, dataNum);// 计算k个分布的协方差矩阵ArrayList<ArrayList<ArrayList<Double>>> covList = new ArrayList<ArrayList<ArrayList<Double>>>();for(int i = 0; i < k; i++) {ArrayList<ArrayList<Double>> dataSetK = new ArrayList<ArrayList<Double>>();for(int j = 0; j < dataNum; j++) {if(type[k] == i) {dataSetK.add(dataSet.get(i));}}covList.set(i, GMMUtil.computeCov(dataSetK, dataDimen, dataSetK.size()));}// 计算每条数据属于每个分布的概率 for(int i = 0; i < dataNum; i++) {for(int j = 0; j < k; j++) {ArrayList<Double> offset = GMMUtil.matrixMinus(dataSet.get(i), pMiu.get(j));ArrayList<ArrayList<Double>> invSigma = covList.get(k);double[] tmp = GMMUtil.matrixSum(GMMUtil.matrixMultiply(GMMUtil.toArray1(offset), GMMUtil.toArray(invSigma)), 2);double coef = Math.pow((2 * Math.PI), -(double)dataDimen / 2d) * Math.sqrt(GMMUtil.computeDet(invSigma, invSigma.size()));px[i][j] = coef * Math.pow(Math.E, -0.5 * tmp[0]);}}return GMMUtil.toList(px);}/** * * @Title: iniParameters * @Description: 初始化参数Parameter* @return Parameter* @throws */public Parameter iniParameters(ArrayList<ArrayList<Double>> dataSet, int dataNum, int k, int dataDimen) {Parameter res = new Parameter();ArrayList<ArrayList<Double>> pMiu = generateCentroids(dataSet, dataNum, k);res.setpMiu(pMiu);// 计算每个样本节点与每个中心节点的距离，以此为据对样本节点进行分类计数，进而初始化k个分布的权值ArrayList<Double> pPi = new ArrayList<Double>();int[] type = getTypes(dataSet, pMiu, k, dataNum);int[] typeNum = new int[k];for(int i = 0; i < dataNum; i++) {typeNum[type[i]]++;}for(int i = 0; i < k; i++) {pPi.add((double)(typeNum[i]) / (double)(dataNum));}res.setpPi(pPi);// 计算k个分布的k个协方差ArrayList<ArrayList<ArrayList<Double>>> pSigma = new ArrayList<ArrayList<ArrayList<Double>>>();for(int i = 0; i < k; i++) {ArrayList<ArrayList<Double>> tmp = new ArrayList<ArrayList<Double>>();for(int j = 0; j < dataNum; j++) {if(type[j] == i) {tmp.add(dataSet.get(i));}}pSigma.add(GMMUtil.computeCov(tmp, dataDimen, dataNum));}res.setpSigma(pSigma);return res;}/** * * @Title: generateCentroids * @Description: 获取随机的k个中心点* @return ArrayList<ArrayList<Double>>* @throws */public ArrayList<ArrayList<Double>> generateCentroids(ArrayList<ArrayList<Double>> data, int dataNum, int k) {ArrayList<ArrayList<Double>> res = null;if(dataNum < k) {return res;} else {res =  new ArrayList<ArrayList<Double>>();List<Integer> random = new ArrayList<Integer>();// 随机产生不重复的k个数while(k > 0) {int index = (int)(Math.random() * dataNum);    if(!random.contains(index)) {    random.add(index);    k--;    res.add(data.get(index));    }}}return res;}/** * * @Title: getTypes * @Description: 返回每条数据的类别* @return int[]* @throws */public int[] getTypes(ArrayList<ArrayList<Double>> dataSet, ArrayList<ArrayList<Double>> pMiu, int k, int dataNum) {int[] type = new int[dataNum];for(int j = 0; j < dataNum; j++) {double minDistance = GMMUtil.computeDistance(dataSet.get(j), pMiu.get(0));type[j] = 0; // 0作为该条数据的类别for(int i = 1; i < k; i++) {if(GMMUtil.computeDistance(dataSet.get(j), pMiu.get(0)) < minDistance) {minDistance = GMMUtil.computeDistance(dataSet.get(j), pMiu.get(0));type[j] = k;}}}return type;}public static void main(String[] args) {ArrayList<Double> pPi = new ArrayList<Double>();System.out.println(pPi.get(0));}}

?一些工具类：

public class GMMUtil {/** * * @Title: computeDistance * @Description: 计算任意两个节点间的距离* @return double* @throws */public static double computeDistance(ArrayList<Double> d1, ArrayList<Double> d2) {double squareSum = 0;for(int i = 0; i < d1.size() - 1; i++) {squareSum += (d1.get(i) - d2.get(i)) * (d1.get(i) - d2.get(i));}return Math.sqrt(squareSum);}/** * * @Title: computeCov * @Description: 计算协方差矩阵* @return ArrayList<ArrayList<Double>>* @throws */public static ArrayList<ArrayList<Double>> computeCov(ArrayList<ArrayList<Double>> dataSet, int dataDimen, int dataNum) {ArrayList<ArrayList<Double>> res = new ArrayList<ArrayList<Double>>();// 计算每一维数据的均值double[] sum = new double[dataDimen];for(ArrayList<Double> data : dataSet) {for(int i = 0; i < dataDimen; i++) {sum[i] += data.get(i);}}for(int i = 0; i < dataDimen; i++) {sum[i] = sum[i] / dataNum;}// 计算任意两组数据的协方差for(int i = 0; i < dataDimen; i++) {ArrayList<Double> tmp = new ArrayList<Double>();for(int j = 0; j < dataDimen; j++) {double cov = 0;for(ArrayList<Double> data : dataSet) {cov += (data.get(i) - sum[i]) * (data.get(j) - sum[j]);}tmp.add(cov);}res.add(tmp);}return res;}/** * * @Title: computeInv * @Description: 计算矩阵的逆矩阵* @return ArrayList<ArrayList<Double>>* @throws */public static double[][] computeInv(ArrayList<ArrayList<Double>> dataSet) {int dataDimen = dataSet.size();double[][] res = new double[dataDimen][dataDimen];// 将list转化为arraydouble[][] a = toArray(dataSet);// 计算伴随矩阵double detA = computeDet(dataSet, dataDimen); // 整个矩阵的行列式 for (int i = 0; i < dataDimen; i++) {              for (int j = 0; j < dataDimen; j++) {                  double num;                  if ((i + j) % 2 == 0) {                      num = computeDet(toList(computeAC(a, i + 1, j + 1)), dataDimen - 1);                  } else {                      num = -computeDet(toList(computeAC(a, i + 1, j + 1)), dataDimen - 1);                  }                  res[j][i] = num / detA;              }          }  return res;}/** * * @Title: computeAC * @Description: 求指定行、列的代数余子式(algebraic complement)* @return double[][]* @throws */public static double[][] computeAC(double[][] dataSet, int r, int c) {int H = dataSet.length;          int V = dataSet[0].length;          double[][] newData = new double[H - 1][V - 1];            for (int i = 0; i < newData.length; i++) {              if (i < r - 1) {                  for (int j = 0; j < newData[i].length; j++) {                      if (j < c - 1) {                          newData[i][j] = dataSet[i][j];                      } else {                          newData[i][j] = dataSet[i][j + 1];                      }                  }              } else {                  for (int j = 0; j < newData[i].length; j++) {                      if (j < c - 1) {                          newData[i][j] = dataSet[i + 1][j];                      } else {                          newData[i][j] = dataSet[i + 1][j + 1];                      }                  }                }          } return newData;}/** * * @Title: computeDet * @Description: 计算行列式* @return double* @throws */public static double computeDet(ArrayList<ArrayList<Double>> dataSet, int dataDimen) {// 将list转化为arraydouble[][] a = toArray(dataSet);if(dataDimen == 2) {return a[0][0] * a[1][1] - a[0][1] * a[1][0];}double res = 0;for(int i = 0; i < dataDimen; i++) {if(i % 2 == 0) {res += a[0][i] * computeDet(toList(computeAC(toArray(dataSet), 1, i + 1)), dataDimen - 1);} else {res += -a[0][i] * computeDet(toList(computeAC(toArray(dataSet), 1, i + 1)), dataDimen - 1);}}return res;}/** * * @Title: toList * @Description: 将array转化成list* @return ArrayList<ArrayList<Double>>* @throws */public static ArrayList<ArrayList<Double>> toList(double[][] a) {ArrayList<ArrayList<Double>> res = new ArrayList<ArrayList<Double>>();for(int i = 0; i < a.length; i++) {ArrayList<Double> tmp = new ArrayList<Double>();for(int j = 0; j < a[i].length; j++) {tmp.add(a[i][j]);}res.add(tmp);}return res;}public static double[][] matrixMultiply(double[][] a, double[][] b) {double[][] res = new double[a.length][b[0].length];for(int i = 0; i < a.length; i++) {for(int j = 0; j < b[0].length; j++) {for(int k = 0; k < a[0].length; k++) {res[i][j] += a[i][k] * b[k][j];}}}return res;}/** * * @Title: dotMatrixMultiply * @Description: 矩阵的点乘，即对应元素相乘* @return double[][]* @throws */public static double[][] dotMatrixMultiply (double[][] a, double[][] b) {double[][] res = new double[a.length][a[0].length];for(int i = 0; i < a.length; i++) {for(int j = 0; j < a[0].length; j++) {res[i][j] = a[i][j] * b[i][j];}}return res;}/** * * @Title: dotMatrixMultiply * @Description: 矩阵的点除，即对应元素相除* @return double[][]* @throws */public static double[][] dotMatrixDivide(double[][] a, double[][] b) {double[][] res = new double[a.length][a[0].length];for(int i = 0; i < a.length; i++) {for(int j = 0; j < a[0].length; j++) {res[i][j] = a[i][j] / b[i][j];}}return res;}/** * * @Title: repmat * @Description: 对应matlab的repmat的函数，对矩阵进行横向或纵向的平铺* @return double[][]* @throws */public static double[][] repmat(double[][] a, int row, int clo) {double[][] res = new double[a.length * row][a[0].length * clo];return null;}/** * * @Title: matrixMinux * @Description: 计算集合只差* @return ArrayList<ArrayList<Double>>* @throws */public static ArrayList<Double> matrixMinus(ArrayList<Double> a1, ArrayList<Double> a2) {ArrayList<Double> res = new ArrayList<Double>();for(int i = 0; i < a1.size(); i++) {res.add(a1.get(i) - a2.get(i));}return res;}/** * * @Title: matrixSum * @Description: 返回矩阵每行之和(mark==2)或每列之和(mark==1)* @return ArrayList<Double>* @throws */public static double[] matrixSum(double[][] a, int mark) {double res[] = new double[a.length];if(mark == 1) { // 计算每列之和，返回行向量res = new double[a[0].length];for(int i = 0; i < a[0].length; i++) {for(int j = 0; j < a.length; j++) {res[i] += a[j][i];}}} else if (mark == 2) { // 计算每行之和， 返回列向量for(int i = 0; i < a.length; i++) {for(int j = 0; j < a[0].length; j++) {res[i] += a[i][j];}}}return res;}public static double[][] toArray(ArrayList<ArrayList<Double>> a) {int dataDimen = a.size();double[][] res = new double[dataDimen][dataDimen];for(int i = 0; i < dataDimen; i++) {for(int j = 0; j < dataDimen; j++) {res[i][j] = a.get(i).get(j);}}return res;}public static double[][] toArray1(ArrayList<Double> a) {int dataDimen = a.size();double[][] res = new double[1][dataDimen];for(int i = 0; i < dataDimen; i++) {res[1][i] = a.get(i);}return res;}/** * * @Title: matrixReverse * @Description: 矩阵专制* @return double[][]* @throws */public static double[][] matrixReverse(double[][] a) {double[][] res = new double[a[0].length][a.length];for(int i = 0; i < a.length; i++) {for(int j = 0; j < a[0].length; j++) {res[j][i] = a[i][j];}}return res;}/** * * @Title: diag * @Description: 向量对角化* @return double[][]* @throws */public static double[][] diag(double[] a) {double[][] res = new double[a.length][a.length];for(int i = 0; i < a.length; i++) {for(int j = 0; j < a.length; j++) {if(i == j) {res[i][j] = a[i];}}}return res;}}

?

?

?