CH9_高斯混合型(GMM)及其Spark实现

Spark自编程实现高斯混合模型(GMM)

1. EM算法简介

EM算法是一种迭代算法，用于含有隐含变量的概率模型参数的极大似然估计，或极大后验概率估计，EM算法由两步组成：E步，求期望；M步，求极大。

2. EM 算法步骤

输入：观测变量数据Y，隐变量数据Z，联合分布P(Y,Z|θ),条件分布P(Z|Y,θ)

输出：模型参数θ

$1.选择参数的初始值{\theta }^{(0)},开始迭代$`

$2.E步：记\theta^{(i)} 为第i次迭代参数 \theta 的估计值,在第i+1次迭代E步，计算：

$$Q(\theta ,{\theta }^{(i)})=E_Z[logP(Y,Z|\theta )|Y,{\theta }^{(i)}]=\sum _{Z}logP(Y,Z|\theta )P(Z|Y,{\theta }^{(i)})$$

$这里P(Z|Y,{\theta }^{(i)}是在给定观测数据Y和当前的参数估计{\theta }^{(i)}下隐变量数据Z的条件概率分布$

$3.M步：求使Q(\theta ,{\theta }^{(i)})极大化的\theta ,确定第i+1次迭代的参数的估计值{\theta }^{(i+1)}$

$${\theta }^{(i+1)}=\underset{\theta }{\underbrace{argmax}}Q(\theta ,{\theta }^{(i)})$$
$4.重复第二步和第三步，直到收敛$

3.EM算法在高斯混合模型(Gaussian Nisture Model)中的应用

高斯混合模型

高斯混合模型是指具有如下形式的概率分布模型：

$$\begin{gathered} P(y|\theta )=\sum _{k=1}^K{\alpha }_k\varphi (y|{\theta }_k) \\ 其中，{\alpha }_k是系数，{\alpha }_k\ge 0,\sum _{k=1}^K{\alpha }_k=1;\varphi (y|{\theta }_k)是高斯分布密度，{\theta }_k=({\mu }_k,{\sigma }_k^2) \end{gathered}$$

$$\varphi (y|{\theta }_k)=\frac{1}{\sqrt{2\pi }{\sigma }_k}-exp(\frac{-(y-{\mu }_k{)}^2}{2{\sigma }_k^2})$$
称为第k个分模型

4. 高斯混合模型参数估计的EM算法

$输入：观测数据y_1,y_2...,y_n，高斯混合模型$

$输出：高斯混合模型参数$

$(1)取参数的初始值开始迭代$

$(2)E步：依据当前模型参数，计算分模型k对观测数据y_i的响应度$

$${\hat{\gamma }}_{jk}=\frac{{\alpha }_k\varphi (y_i|{\theta }_k)}{{\sum }_{k=1}^K{\alpha }_k\varphi (y_i|{\theta }_k)}j=1,2,...,N;k=1,2,...,K$$

$(3)M步：计算新一轮迭代的模型参数$

$${\hat{\mu }}_k=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}{y}_j}{\sum {j=1}^N\hat{\gamma }jk},k=1,2,...,k$$

$${\hat{\sigma }}_k^2=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}(y_i-{\mu }_k{)}^2}{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}},k=1,2...,K$$

$${\hat{\alpha }}_k=\frac{{\sum }_{j=1}^N{\hat{\gamma }}_{jk}}{N},k=1,2...,K$$

$(4)重复第二步和第三步，知道收敛或达到最大迭代次数$

---

Spark自编程实现GMM

package CH9_EM

import org.apache.spark.sql.{DataFrame, Row, SparkSession}
import org.apache.spark.sql.functions._
import org.apache.spark.ml.linalg.{DenseMatrix, DenseVector, Vector}

import scala.beans.BeanProperty
import breeze.linalg.{
  *,
  diag,
  DenseMatrix => denseMatrix,
  DenseVector => densevector
}
import org.apache.spark.ml.feature.VectorAssembler
import org.apache.spark.ml.stat.distribution.MultivariateGaussian
import org.apache.spark.ml.util.Identifiable
import org.apache.spark.rdd.RDD
import org.apache.spark.sql.types.{IntegerType, StructType}

/**
  * Created by WZZC on 2020/5/17
  **/
case class GMMModel(data: DataFrame) {

  @BeanProperty var maxIter: Int = 40
  @BeanProperty var tol: Double = 1e-5
  @BeanProperty var k: Int = _
  @BeanProperty var fts: Array[String] = data.columns

  //  private val spark: SparkSession = data.sparkSession

  private val numFeatures = fts.length
  private val dataSize = data.count()

  // 初始化 alpha、 mu和sigma向量

  var alpha: Array[Double] = _

  var mu: Array[densevector[Double]] = _

  var cov: Array[denseMatrix[Double]] = _

  private def matrixArr: densevector[Double] =
    densevector(densevector.rand(numFeatures).toArray)
      .asInstanceOf[densevector[Double]]

  // 结果表的Schema信息
  val schemaOfResult: StructType = data.schema
    .add("clusterId", IntegerType) //增加一列表示类id的字段

  private var clusterDf: DataFrame = _
//  var clusterDf  = spark.createDataFrame(spark.sparkContext.emptyRDD[Row], schemaOfResult)

  private val ftsName: String = Identifiable.randomUID("GMMModel")

  /**
    * 数据特征转换
    *
    * @param dataFrame
    * @return
    */
  def dataTransForm(dataFrame: DataFrame) = {
    new VectorAssembler()
      .setInputCols(fts)
      .setOutputCol(ftsName)
      .transform(dataFrame)
  }

  /**
    * 计算多元高斯分布的概率密度函数
    *
    * @param vec 待计算样本
    * @param mus 均值向量
    * @param covs 协方差矩阵
    * @return
    */
  private def pdf(vec: densevector[Double],
                  mus: densevector[Double],
                  covs: denseMatrix[Double]) = {

    val muvec: DenseVector = new DenseVector(mus.data)
    val yvec: DenseVector = new DenseVector(vec.data)
    val covMatrix =
      new DenseMatrix(covs.rows, covs.cols, covs.data, covs.isTranspose)

    val mgaussian: MultivariateGaussian =
      new MultivariateGaussian(muvec, covMatrix)
    //    val mgaussian: MultivariateGaussian = MultivariateGaussian(mus, covs)
    //    使用breeze包的多元高斯分布会出现以下错误
    //    breeze.linalg.NotConvergedException
    mgaussian.pdf(yvec)

  }

  /**
    *
    * @return
    */
  private def gamakudf(mus: Array[densevector[Double]],
                       covs: Array[denseMatrix[Double]],
                       alphas: Array[Double]) =
    udf((vec: Vector) => {

      val tuples: Array[((densevector[Double], denseMatrix[Double]), Double)] =
        mus.zip(covs).zip(alphas)

      val featureVec = new densevector(vec.toArray)

      val gammak: Array[Double] = tuples.map(tp => {
        pdf(featureVec, tp._1._1, tp._1._2) * tp._2
      })

      val s = gammak.sum

      new DenseVector(gammak.map(x => (x / s).formatted("%.4f").toDouble))
 
    })

 

  private def pdfudf(mus: Array[densevector[Double]],
                     covs: Array[denseMatrix[Double]]) =
    udf((vec: Vector) => {

      val tuples = mus.zip(covs)

      val featureVec = new densevector(vec.toArray)

      val gammak: Array[Double] = tuples.map(tp => {
        pdf(featureVec, tp._1, tp._2)
      })

      new DenseVector(gammak)

    })

  /**
    *
    * @param matrix
    * @param side
    */
  private def matrixToVectors(matrix: denseMatrix[Double], side: String) = {

    val rows: Int = matrix.rows
    val cols: Int = matrix.cols

    side match {
      case "row" => {
        (0 until rows).toArray
          .map(i => {
            densevector((0 until cols).toArray.map(j => {
              matrix.valueAt(i, j)
            }))
          })
      }

      case "col" => {
        (0 until cols).toArray
          .map(i => {
            densevector((0 until rows).toArray.map(j => {
              matrix.valueAt(j, i)
            }))
          })
      }
    }

  }

  def fit = {

    var alphas: Array[Double] = new Array[Double](k).map(_ => 1.0 / k)

    var mus: Array[densevector[Double]] =
      Array.fill(k)(densevector.rand(numFeatures))

    var covs: Array[denseMatrix[Double]] =
      Array.fill(k)(diag(matrixArr))

    var i = 0
//    var cost = 0.0
//    var convergence = false //判断收敛，即代价函数变化小于阈值tol

    while (i < maxIter) {

      //  E步:根据当前模型参数，计算分模型k对观测数据yi的响应度
      val edf: DataFrame = dataTransForm(data)
        .withColumn("gammajk", gamakudf(mus, covs, alphas)(col(ftsName)))

      val gammajk: RDD[densevector[Double]] = edf
        .select("gammajk")
        .rdd
        .map(x => densevector(x.getAs[Vector](0).toArray))

      //  M步:计算新一轮迭代的模型参数
      val gammaMatrix: denseMatrix[Double] = gammajk
        .map(_.toDenseMatrix)
        .reduce((v1, v2) => denseMatrix.vertcat(v1, v2))
        .t

      val sumgammak: Array[Double] = gammajk
        .map(_.toArray)
        .reduce((a1, a2) => a1.zip(a2).map(x => x._1 + x._2))

      val dataRdd = dataTransForm(data)
        .select(ftsName)
        .rdd
        .persist()

      val gammaVectors: Array[densevector[Double]] =
        matrixToVectors(gammaMatrix, "row")

      val dfmatrix: denseMatrix[Double] = dataRdd
        .map(x => densevector(x.getAs[Vector](0).toArray))
        .map(_.toDenseMatrix)
        .reduce((v1, v2) => denseMatrix.vertcat(v1, v2))

      //
      // 更新sigmas
      covs = mus
        .zip(gammaVectors)
        .map(tp => {

          val ymusMatrix: denseMatrix[Double] = dataRdd
            .map(x => {
              val ym = densevector(x.getAs[Vector](0).toArray)
              ym - tp._1
            })
            .map(_.toDenseMatrix)
            .reduce((v1, v2) => denseMatrix.vertcat(v1, v2))

          (ymusMatrix(::, *) *:* tp._2).t * ymusMatrix

        })
        .zip(sumgammak)
        .map(tp => {
          tp._1.map(x => (x / tp._2).formatted("%.4f").toDouble)
        })

      
      // 更新 mus
      mus = matrixToVectors(gammaMatrix, "row")
        .map(vec => {
          (vec.toDenseMatrix * dfmatrix).toDenseVector
        })
        .zip(sumgammak)
        .map(tp => {
          tp._1.map(_ / tp._2)
        })

      // 更新 alphas
      alphas = sumgammak.map(_ / dataSize).map(_.formatted("%.4f").toDouble)

      i += 1
      if (i >= maxIter) clusterDf = edf
    }

    alpha = alphas
    mu = mus
    cov = covs

  }

  //TODO   predict
  def predict = {

    val predictUdf = udf((probability: Vector) => {
      val array = probability.toArray
      array.indexOf(array.max)
    })

    clusterDf
      .withColumnRenamed("gammajk", "probability")
      .withColumn("prediction", predictUdf(col("probability")))
      .drop(ftsName)

  }

}

 

模型测试

bject GMMModelRunner {

  def main(args: Array[String]): Unit = {

    val spark = SparkSession
      .builder()
      .appName(s"${this.getClass.getSimpleName}")
      .master("local[*]")
      .config("spark.serializer", "org.apache.spark.serializer.KryoSerializer")
      .config("spark.shuffle.consolidateFiles", "true")
      .config("spark.io.compression.codec", "snappy")
      .getOrCreate()

    val dataset: DataFrame = spark.read
      .option("header", true)
      .option("inferSchema", true)
      .csv("F:\\DataSource\\gmm.txt")

    val fts = Array("f1", "f2", "f3")

    val model: GMMModel = new GMMModel(dataset)
    model.setFts(fts)
    model.setK(2)
    model.setMaxIter(5)

    model.fit

    val alphas = model.alpha
    val mus = model.mu
    val covs = model.cov

    println("================")
    alphas.foreach(a => println("alpha: " + a))
    println("================")
    covs.foreach(println)
    println("================")
    mus.foreach(println)

    model.predict.show(false)
    spark.stop()

  }

}

参考资料

《统计学习方法》