【Scala统计学系列】scala Spark breeze中Gamma伽马分布 ChiSquared卡方分布 用途方法示例源码分析

Spark breeze中Gamma伽马分布 ChiSquared卡方分布 Gaussian高斯分布 用途方法示例源码分析

用途

伽马分布在统计学和概率论中有广泛的应用。在Spark中,Gamma分布主要用于以下几个方面:

  1. 数据生成和模拟:通过使用Gamma分布生成符合特定形状和尺度参数的随机数,可以模拟各种类型的数据。这对于测试和验证数据处理流程以及算法的性能非常有用。

  2. 统计分析:Gamma分布在统计学中具有广泛的应用。使用Gamma分布可以对数据进行建模和估计,例如拟合观测数据、计算概率密度函数计算累积分布函数等。Spark中的Gamma分布类提供了计算这些统计量的方法。

  3. 机器学习和数据挖掘:Gamma分布在一些机器学习和数据挖掘算法中起到重要作用,例如在概率图模型中表示变量之间的依赖关系,或者作为损失函数的一部分来优化模型参数。

  4. 随机采样和模型评估:在某些情况下,需要从给定的Gamma分布中采样数据,用于构建训练集、验证集或测试集。此外,使用Gamma分布生成的随机样本还可以用于评估模型的性能和鲁棒性。

总而言之,Spark中的Gamma分布类提供了方便的方法来处理和操作Gamma分布相关的数据和问题,无论是进行统计分析、模拟数据还是应用于机器学习和数据挖掘任务。

特例(卡方分布)

伽马分布与卡方分布的关系:

当卡方分布的自由度为正整数n时,它可以表示为一个形状参数为n/2、尺度参数为2的伽马分布。即:

卡方分布(X ~ χ²(n)) = 伽马分布(X ~ Γ(n/2, 2))

这意味着从伽马分布中采样得到的随机变量经过适当的标准化处理后可以得到卡方分布。具体而言,如果从伽马分布中采样得到的随机变量X,将其除以自由度n,即 X/n,则得到的结果符合自由度为n的卡方分布。注意:当自由度参数为整数时,卡方分布是伽马分布的特例。

在统计学中,卡方分布是一种特殊的伽马分布,它描述了自由度为正整数独立标准正态随机变量的平方和。如果将一个服从自由度为n的卡方分布的随机变量X进行标准化处理,即将其除以n,那么得到的结果服从形状参数为n/2、尺度参数为2的伽马分布。

这个关系非常有用,因为卡方分布在统计推断和假设检验中经常出现。通过理解伽马分布和卡方分布之间的关系,可以更好地理解和应用卡方分布,并将其与其他概率分布进行比较和联系。

示例1-gamma

object GammaTest extends App{
  import breeze.stats.distributions.Gamma

  // 创建一个Gamma分布对象
  val gammaDist = Gamma(shape = 2.0, scale = 3.0)

  // 计算概率密度函数值
  val pdfValue = gammaDist.pdf(4.0)

  // 从Gamma分布中抽取值
  val drawValue = gammaDist.draw()

  // 计算均值
  val meanValue = gammaDist.mean

  // 计算方差
  val varianceValue = gammaDist.variance

  // 计算众数
  val modeValue = gammaDist.mode

  // 计算熵
  val entropyValue = gammaDist.entropy

  // 计算累积分布函数值
  val cdfValue = gammaDist.cdf(3.0)
}

请注意,上述代码示例仅为演示如何使用breeze库中的Gamma类创建和操作伽马分布。您需要确保已经正确导入breeze.stats.distributions.Gamma类,并根据自己的需求调整参数和方法的使用。

示例2-chi

在Spark中,可以使用breeze.stats.distributions.Gammabreeze.stats.distributions.ChiSquared类来生成伽马分布和卡方分布,并计算其概率密度函数(PDF)的值。下面是一个简单的代码示例:

object chi2GammaTest extends App{
  import breeze.stats.distributions.{Gamma, ChiSquared}

  // 创建伽马分布对象
  val gammaDist = Gamma(shape = 2.0, scale = 3.0)

  // 创建卡方分布对象
  val chiSquareDist = ChiSquared(4)

  // 计算伽马分布和卡方分布的概率密度函数值
  val x = 2.5
  val gammaPdf = gammaDist.pdf(x)
  val chiSquarePdf = chiSquareDist.pdf(x)

  println(s"伽马分布的 PDF 值:$gammaPdf")
  println(s"卡方分布的 PDF 值:$chiSquarePdf")
//  伽马分布的 PDF 值:0.12072172458529948
//  卡方分布的 PDF 值:0.17906549803761884
}

示例3-gauss

package org.example.scala

import breeze.linalg.DenseVector
import breeze.stats.distributions.Gaussian

object GaussianDemo {
  def main(args: Array[String]): Unit = {
    val mean = 0.0 // 均值
    val stddev = 1.0 // 标准差

    val gaussian = new Gaussian(mean, stddev) // 创建高斯分布对象

    val x = DenseVector(-2.0, -1.0, 0.0, 1.0, 2.0) // 输入数据

    // 计算每个输入数据对应的高斯函数值
    val probabilities = x.map(gaussian.pdf)

    println("输入数据:")
    x.foreach(println)

    println("\n高斯函数值:")
    probabilities.foreach(println)
  }
}
//输入数据:
//-2.0
//-1.0
//0.0
//1.0
//2.0
//
//高斯函数值:
//0.053990966513188056
//0.24197072451914337
//0.3989422804014327
//0.24197072451914337
//0.053990966513188056

方法

pdf 计算概率密度函数值。

override def pdf(x: Double): Double

参数:

  • x:变量x

返回值:概率密度值

logDraw从伽马分布中抽取对数值。

def logDraw(): Double

返回值:对数值

draw 从伽马分布中抽取值。

def draw(): Double

返回值:值

mean 计算均值。

def mean: Double

返回值:均值

variance 计算方差。

def variance: Double

返回值:方差

mode计算众数。

def mode: Double

返回值:众数

entropy计算熵。

def entropy: Double

返回值:熵

probability计算概率。

override def probability(x: Double, y: Double): Double

参数:

  • x:变量x
  • y:变量y

返回值:概率值

inverseCdf计算反函数。

override def inverseCdf(p: Double): Double

参数:

  • p:概率值

返回值:反函数值

cdf计算累积分布函数值。

override def cdf(x: Double): Double

参数:

  • x:变量x

返回值:累积分布函数值

伴生对象

object Gamma extends ExponentialFamily[Gamma, Double] with ContinuousDistributionUFuncProvider[Double, Gamma]

伽马分布的伴生对象,提供了一些额外的方法和函数。

SufficientStatistic 伽马分布的充分统计量。

case class SufficientStatistic(n: Double, meanOfLogs: Double, mean: Double) extends BaseSuffStat[SufficientStatistic]

emptySufficientStatistic创建一个空的充分统计量。

def emptySufficientStatistic: SufficientStatistic

返回值:空的充分统计量

sufficientStatisticFor根据给定的数据计算充分统计量。

def sufficientStatisticFor(t: Double): SufficientStatistic

参数:

  • t:数据

返回值:充分统计量

mle最大似然估计。

def mle(ss: SufficientStatistic): (Double, Double)

参数:

  • ss:数据的充分统计量

返回值:形状参数和尺度参数的元组

distribution根据给定的参数创建一个伽马分布对象。

def distribution(p: Parameter): Gamma

参数:

  • p:参数

返回值:伽马分布对象

likelihoodFunction计算似然函数。

def likelihoodFunction(stats: SufficientStatistic): DiffFunction[Parameter]

参数:

  • stats:数据的充分统计量

返回值:似然函数

伽马中文源码


/**
 * Represents a Gamma distribution.
 * E[X] = shape * scale
 *
 * @author dlwh
 */
case class Gamma(shape : Double, scale : Double)(implicit rand: RandBasis = Rand)
    extends ContinuousDistr[Double] with Moments[Double, Double] with HasCdf with HasInverseCdf {
  if(shape <= 0.0 || scale <= 0.0)
    throw new IllegalArgumentException("Shape and scale must be positive")

  override def pdf(x:Double) = if (x > 0) {
    math.exp(logPdf(x))
  } else {
    if (shape > 1.0) {
      0.0
    } else if (shape == 1.0) {
      normalizer
    } else {
      Double.PositiveInfinity
    }
  }

  lazy val logNormalizer: Double = lgamma(shape) + shape * log(scale)

  override def unnormalizedLogPdf(x : Double) = (shape - 1) * log(x) - x/scale

  override def toString = "Gamma(" + shape + "," + scale + ")"

  def logDraw() = if (shape < 1) {
    // adapted from numpy distributions.c which is Copyright 2005 Robert Kern (robert.kern@gmail.com) under BSD
    @tailrec
    def rec: Double = {
      val u = rand.uniform.draw()
      val v = -math.log(rand.uniform.draw())
      val logU = log(u)
      if ( logU <= math.log1p(-shape)) {
        val logV = log(v)
        val logX = logU / shape
        if (logX <= logV)  logX
        else rec
      } else {
        val y = -log((1-u)/shape)
        val logX = math.log(1.0 - shape + shape*y)/shape
        if (logX <= math.log(v + y)) logX
        else rec
      }
    }
    rec + math.log(scale)
  } else math.log(draw)

  def draw() = {
    if(shape == 1.0) {
      scale * -math.log(rand.uniform.draw())
    } else if (shape < 1.0) {
      // from numpy distributions.c which is Copyright 2005 Robert Kern (robert.kern@gmail.com) under BSD
      @tailrec
      def rec: Double = {
        val u = rand.uniform.draw()
        val v = -math.log(rand.uniform.draw())
        if (u <= 1.0 - shape) {
          val x = pow(u, 1.0 / shape)
          if (x <= v)  x
          else rec
        } else {
          val y = -log((1-u)/shape)
          val x = pow(1.0 - shape + shape*y, 1.0 / shape)
          if (x <= (v + y)) x
          else rec
        }
      }

      scale * rec
//      val c = 1.0 + shape/scala.math.E
//      var d = c * rand.uniform.draw()
//      var ok = false
//      var x = 0.0
//      while(!ok) {
//        if (d >= 1.0) {
//          x = -log((c - d) / shape)
//          if (-math.log(rand.uniform.draw()) >= (1.0 - shape) * log(x)) {
//            x = (scale * x)
//            ok = true
//          }
//        } else {
//          x = math.pow(d, 1.0/shape)
//          if (-math.log(rand.uniform.draw()) >= (1.0 - shape) * log(x)) {
//            x = scale * x
//            ok = true
//          }
//        }
//        d = c * rand.uniform.draw()
//      }
//      x
    } else {
      // from numpy distributions.c which is Copyright 2005 Robert Kern (robert.kern@gmail.com) under BSD
      val d = shape-1.0/3.0
      val c = 1.0 / math.sqrt(9.0* d)
      var r = 0.0
      var ok = false
      while (!ok) {
        var v = 0.0
        var x = 0.0
        do {
          x = rand.generator.nextGaussian
          v = 1.0 + c * x
        } while(v <= 0)

        v = v*v*v
        val x2 = x * x
        val u = rand.uniform.draw()
        if (  u < 1.0 - 0.0331 * (x2 * x2)
          || log(u) < 0.5*x2 + d* (1.0 - v+log(v))) {
          r = (scale*d*v)
          ok = true
        }
      }
      r
    }
  }

  def mean = shape * scale
  def variance = mean * scale
  def mode = { require(shape >= 1); mean - scale}
  def entropy = logNormalizer - (shape - 1) * digamma(shape) + shape

  override def probability(x: Double, y: Double): Double = {
    new GammaDistribution(shape, scale).probability(x, y)
  }

  override def inverseCdf(p: Double): Double = {
//    gammp(this.shape, p / this.scale);
    new GammaDistribution(shape, scale).inverseCumulativeProbability(p)
  }

  override def cdf(x: Double): Double = {
    new GammaDistribution(shape, scale).cumulativeProbability(x)
  }
}

object Gamma extends ExponentialFamily[Gamma,Double] with ContinuousDistributionUFuncProvider[Double,Gamma] {
  type Parameter = (Double,Double)
  import breeze.stats.distributions.{SufficientStatistic => BaseSuffStat}
  case class SufficientStatistic(n: Double, meanOfLogs: Double, mean: Double) extends BaseSuffStat[SufficientStatistic] {
    def *(weight: Double) = SufficientStatistic(n*weight, meanOfLogs, mean)
    def +(t: SufficientStatistic) = {
      val delta = t.mean - mean
      val newMean = mean + delta * (t.n /(t.n + n))
      val logDelta = t.meanOfLogs - meanOfLogs
      val newMeanLogs = meanOfLogs + logDelta * (t.n /(t.n + n))
      SufficientStatistic(t.n+n, newMeanLogs, newMean)
    }
  }

  def emptySufficientStatistic = SufficientStatistic(0,0,0)

  def sufficientStatisticFor(t: Double) = SufficientStatistic(1,math.log(t),t)

  // change from Timothy Hunter. Thanks!
  def mle(ss: SufficientStatistic) = {
    val s = math.log( ss.mean ) - ss.meanOfLogs
    assert(s > 0 , s) // check concavity
    val k_approx = approx_k(s)
    assert(k_approx > 0 , k_approx)
    val k = Nwt_Rph_iter_for_k(k_approx, s)
    val theta = ss.mean / (k)
    (k, theta)
  }
  /*
   * s = log( x_hat) - log_x_hat
   */
  def approx_k(s:Double) : Double = {
    // correct within 1.5%
    (3 - s + math.sqrt( math.pow((s-3),2) + 24*s )) / (12 * s)
  }

  def Nwt_Rph_iter_for_k(k:Double, s:Double ): Double = {
    /*
     * For a more precise estimate, use Newton-Raphson updates
     */
    val k_new = k - (math.log(k) - digamma(k) - s)/( 1.0/k - trigamma(k) )
    if (math.abs(k - k_new)/math.abs(k_new) < 1.0e-4)
      k_new
    else
      Nwt_Rph_iter_for_k(k_new,s)
  }


  def distribution(p: Parameter) = new Gamma(p._1,p._2)

  def likelihoodFunction(stats: SufficientStatistic) = new DiffFunction[(Double,Double)] {
    val SufficientStatistic(n,meanOfLogs,mean) = stats
    def calculate(x: (Double,Double)) = {
      val (a,b) = x
      val obj = -n * ((a - 1) * meanOfLogs - lgamma(a) - a * log(b)- mean / b)
      val gradA = -n * (meanOfLogs - digamma(a)  - log(b))
      val gradB = -n * (-a/b + mean / b / b)
      (obj,(gradA,gradB))
    }
  }
}

高斯中文源码

/**
 * 表示单个实数变量的高斯分布(正态分布)。
 *
 * @param mu 均值
 * @param sigma 标准差
 * @param rand 随机数生成器
 */
case class Gaussian(mu: Double, sigma: Double)(implicit rand: RandBasis = Rand)
    extends ContinuousDistr[Double] with Moments[Double, Double] with HasCdf with HasInverseCdf {
  // 使用随机数生成器生成符合指定均值和标准差的高斯分布随机数
  private val inner = rand.gaussian(mu, sigma)

  /**
   * 从高斯分布中抽取一个随机数。
   *
   * @return 随机数
   */
  def draw(): Double = inner.get()

  override def toString(): String = "Gaussian(" + mu + ", " + sigma + ")"

  /**
   * 计算高斯分布变量在区间[x, y]内的概率。
   * 这个概率被计算为 P[Z < y] - P[Z < x],其中Z是高斯随机变量。
   *
   * @param x 区间的下限
   * @param y 区间的上限
   * @return 高斯随机变量Z落在区间[x, y]内的概率
   */
  override def probability(x: Double, y: Double): Double = {
    require(x <= y, "Undefined probability: lower-end of the interval should be smaller than its upper-end")
    cdf(y) - cdf(x)
  }

  /**
   * 计算高斯分布的反函数,即给定概率值p,计算使得累积分布函数(CDF)等于p的对应变量x。
   *
   * @param p 概率值在[0,1]范围内
   * @return CDF等于p时对应的变量x
   */
  def inverseCdf(p: Double): Double = {
    require(p >= 0)
    require(p <= 1)

    mu + sigma * Gaussian.sqrt2 * erfinv(2 * p - 1)
  }

  /**
   * 计算高斯分布的累积分布函数(CDF)。
   *
   * @param x 变量x
   * @return 累积分布函数值
   */
  def cdf(x: Double): Double = .5 * (1 + erf((x - mu) / (Gaussian.sqrt2 * sigma)))

  override def unnormalizedLogPdf(t: Double): Double = {
    val d = (t - mu) / sigma
    -d * d / 2.0
  }

  override lazy val normalizer: Double = 1.0 / sqrt(2 * Pi) / sigma
  lazy val logNormalizer: Double = log(sqrt(2 * Pi)) + log(sigma)

卡方中文源码

/**
 * 表示自由度为k的卡方分布。
 *
 * @param k 自由度参数
 * @param rand 随机数生成器
 */
case class ChiSquared(k: Double)(implicit rand: RandBasis = Rand)
    extends ContinuousDistr[Double] with Moments[Double, Double] with HasCdf with HasInverseCdf {
  // 内部使用伽马分布来实现卡方分布
  private val innerGamma = Gamma(k / 2, 2)

  /**
   * 从卡方分布中抽取一个随机数。
   *
   * @return 随机数
   */
  def draw(): Double = innerGamma.draw()

  /**
   * 计算卡方分布的概率密度函数(PDF)。
   *
   * @param x 变量x
   * @return 概率密度函数值
   */
  override def pdf(x: Double): Double = if (x > 0.0) {
    math.exp(logPdf(x))
  } else if (x == 0.0) {
    if (k > 2.0) {
      0.0
    } else if (k == 2.0) {
      0.5
    } else {
      Double.PositiveInfinity
    }
  } else {
    throw new IllegalArgumentException("Domain of ChiSquared.pdf is [0,Infinity), you tried to apply to " + x)
  }

  /**
   * 计算未归一化的对数概率密度函数值。
   *
   * @param x 变量x
   * @return 未归一化的对数概率密度函数值
   */
  def unnormalizedLogPdf(x: Double): Double = innerGamma.unnormalizedLogPdf(x)

  /**
   * 获取卡方分布的归一化常数的对数值。
   *
   * @return 归一化常数的对数值
   */
  lazy val logNormalizer: Double = innerGamma.logNormalizer

  /**
   * 计算卡方分布的均值。
   *
   * @return 均值
   */
  def mean: Double = innerGamma.mean

  /**
   * 计算卡方分布的方差。
   *
   * @return 方差
   */
  def variance: Double = innerGamma.variance

  /**
   * 计算卡方分布的众数。
   *
   * @return 众数
   */
  def mode: Double = innerGamma.mode

  /**
   * 计算卡方分布的熵。
   *
   * @return 熵
   */
  def entropy: Double = innerGamma.entropy

  override def toString: String = ScalaRunTime._toString(this)

  /**
   * 计算卡方分布的概率。
   *
   * @param x 变量x
   * @param y 变量y
   * @return 概率值
   */
  override def probability(x: Double, y: Double): Double = {
    innerGamma.probability(x, y)
  }

  /**
   * 计算卡方分布的反函数。
   *
   * @param p 概率值
   * @return 反函数值
   */
  override def inverseCdf(p: Double): Double = {
    innerGamma.inverseCdf(p)
  }

  /**
   * 计算卡方分布的累积分布函数值。
   *
   * @param x 变量x
   * @return 累积分布函数值
   */
  override def cdf(x: Double): Double = {
    innerGamma.cdf(x)
  }
}

object ChiSquared extends ExponentialFamily[ChiSquared, Double] with ContinuousDistributionUFuncProvider[Double, ChiSquared] {
  type Parameter = Double
  type SufficientStatistic = Gamma.SufficientStatistic

  /**
   * 创建一个空的充分统计量。
   *
   * @return 空的充分统计量
   */
  def emptySufficientStatistic: ChiSquared.SufficientStatistic = Gamma.emptySufficientStatistic

  /**
   * 根据给定的数据计算充分统计量。
   *
   * @param t 数据
   * @return 充分统计量
   */
  def sufficientStatisticFor(t: Double): ChiSquared.SufficientStatistic = Gamma.sufficientStatisticFor(t)

  /**
   * 计算卡方分布的似然函数。
   *
   * @param stats 数据的充分统计量
   * @return 似然函数
   */
  def likelihoodFunction(stats: ChiSquared.SufficientStatistic): DiffFunction[ChiSquared.Parameter] = {
    val inner = Gamma.likelihoodFunction(stats)
    new DiffFunction[ChiSquared.Parameter] {
      def calculate(x: ChiSquared.Parameter): (Double, ChiSquared.Parameter) = {
        val (obj, ggrad) = inner.calculate((x / 2.0, 2.0))
        obj -> ggrad._1
      }
    }
  }

  /**
   * 使用最大似然估计(MLE)计算卡方分布的参数。
   *
   * @param ss 数据的充分统计量
   * @return 参数值
   */
  def mle(ss: ChiSquared.SufficientStatistic): ChiSquared.Parameter = {
    ss.mean
  }

  /**
   * 创建一个卡方分布对象。
   *
   * @param p 参数
   * @return 卡方分布对象
   */
  def distribution(p: ChiSquared.Parameter): ChiSquared = ChiSquared(p)
}

上述代码是在Spark中定义的ChiSquared类和伴生对象,并添加了中文注释。ChiSquared类表示自由度为k的卡方分布,提供了计算概率密度函数、均值、方差等统计量的方法。伴生对象实现了最大似然估计(MLE)和其他相关功能。

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

BigDataMLApplication

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值