离散概率类套件

Joe Nellis

4.93/5 (19投票s)

2003年3月27日

9分钟阅读

73261

2215

用于计算离散概率的五个类。

摘要

本文介绍五种最常见的离散概率分布，这些分布可用于计算机模拟、游戏、人工智能决策制定和环境建模。涵盖的分布包括二项分布、几何分布、负二项分布、超几何分布和泊松概率分布，以及一个抽象的概率类。这些类与其他常规的概率计算数学方法不同之处在于，它们采用了几项实现技巧，扩展了允许输入的范围。文章还对蒙特卡洛方法与概率计算进行了比较。

概述

如果你是计算机科学或数学专业的学生，那么你很可能会上概率与统计课程。概率是一个庞大的研究领域，但基本概念任何人都可以学习。以下五个概率类很可能是你在学习中最先遇到的，它们属于“离散”概率的范畴，这意味着它们模拟的事件要么发生，要么不发生。我将简要介绍每种概率及其可能的应用场景。离散概率涉及大量阶乘，或者通常需要将非常大的数与非常小的数相乘。文中讨论了这些问题以及解决方法，并介绍了许多程序员用来通过随机数生成概率的常用替代方法。

概率基础

概率是指事件发生的可能性。概率始终在 0.0 到 1.0 之间（包括 0.0 和 1.0）。概率总是可以定义为“在一个实验中事件发生的次数”除以“实验所有可能结果的总次数”。掷骰子是除抛硬币之外最常见的例子。掷出任何一个数字的概率是 1 除以骰子的面数。对于普通骰子，通常是 1/6。抛硬币出现正面的概率是 1/2。这些都是显而易见的观察结果，但它们是更深层次问题的基础，例如连续掷三次骰子出现相同数字的概率是多少。这并不那么明显。我们将要看到的这五种概率分布都依赖于对离散事件的先验知识，例如掷骰子、抛硬币、击球、射击目标、设备故障、彩票或钻井取油等。如果我们知道离散事件发生的可能性，我们就可以开始询问在特定情况下事件如何发展的问题。

二项概率

二项概率模型用于计算在总试验次数 (N) 中，成功试验的次数 Y（也称为“随机变量”）的概率。每次试验成功的几率（或失败的几率）都完全相同。例如，掷骰子，掷出 1 的概率是 1/6。如果我们说成功（用 p 表示）p=1/6，那么掷出 1 的概率就是 1.0 - p，用 q 表示。符号 p 表示成功，符号 q 表示失败。假设我们有五个骰子，或者我们掷同一个骰子 5 次，并且我们想知道在 5 次掷骰子中恰好出现 3 次 1 的概率。这里的 3 是成功掷出 1 的次数。二项概率的完整公式是

P(Y) = N! / ((N-Y)! * Y!)  * pY * q (N-Y)

在这个例子中，N = 5，Y = 3，p=1/6，q=5/6。该概率表示为：P(Y=3)。如果我们想知道在五次掷骰子中出现“至少三个”或“三个或更少”的 1 的概率，我们将说 P(Y<=3)。这称为累积概率，简单地是 P(Y=0)+P(Y=1)+P(Y=2)+P(Y=3) 的总和。抽象类 Probability 包含一个名为 m_RVT 的成员变量，意为“随机变量类型”，它是一个枚举，用于指示 P(Y ? y) 的符号。GetResult 方法会查看此枚举，并在派生类中为累积概率中的每个 P(Y=y) 调用 ComputeResult 函数。

double Probability::GetResult()  throw (ProbabilityException)
{
  if(m_probability != UNDEFINED )
      return m_probability;
  try{
      m_probability = 0.0;
    int i = 0;
    switch(m_RVT)
    {
    case EQUAL:
        m_probability = ComputeResult();
        break;
    case LESSTHAN:
        for(SetRV(i); i < m_RV; SetRV(++i))
            m_probability += ComputeResult();
        break;
    case GREATERTHANOREQUAL:
        for(SetRV(i),m_probability = 1.0; i < m_RV; SetRV(++i))
            m_probability -= ComputeResult();
        break;
    case GREATERTHAN:
        for(SetRV(i),m_probability = 1.0; i <= m_RV; SetRV(++i))
            m_probability -= ComputeResult();
        break;
    case LESSTHANOREQUAL:
        for(SetRV(i); i <= m_RV; SetRV(++i))
            m_probability += ComputeResult();
        break;
    case NOTEQUAL:
        m_probability = 1.0 - ComputeResult();
        return m_probability;
    default:
        throw ProbabilityException(
            "Probability error- Random Variable type is unset");
    }
  }
  catch(ProbabilityException pe)
  {
    m_probability = UNDEFINED;
    SetRV(m_RV);
    throw pe;
  }
  SetRV(m_RV);
  return m_probability;
}

抽象类 Probability 的定义很小且相当直接

class Probability  
{
public:
  enum RandomVariableType { EQUAL, LESSTHAN, GREATERTHAN, 
                  LESSTHANOREQUAL, GREATERTHANOREQUAL, NOTEQUAL };

  /* Constructor Parameters:
    rv - The Random Variable value.
    rvt - Random Variable type, whether it is cumulative 
          and which way it is cumulative.
    prob - The probability of the event.  This can be set before 
           hand if it is known and the GetResult function will return it. 
  */
  Probability( int rv=0, RandomVariableType rvt=EQUAL, double prob=UNDEFINED)
    :m_RV(rv), m_RVT(rvt), m_probability(prob){};
  virtual ~Probability(){};

  class ProbabilityException  
  {
  public:
    ProbabilityException(const char* whatString):What(whatString){}
    inline const char* what() const { return What; }
  protected:
    const char* What;
  };

  virtual double GetResult() throw (ProbabilityException) ;
  virtual double GetExpectedValue() const = 0;
  virtual double GetVariance() const = 0;
protected:
  virtual double ComputeResult() = 0;
  virtual void   SetRV(int Y)    = 0;
  RandomVariableType m_RVT;
  double m_probability;
  int m_RV;
};

每个专门的概率类构造函数都会像这样调用此基类构造函数

BinomialProbability::BinomialProbability(int N , int Y, double p,
                                         Probability::RandomVariableType rvt)
        :m_trials(N),
         m_successes(Y),
         m_chance_of_success(p),
        Probability(Y, rvt)
{
    assert(Y<=N && Y >=0);
    assert(p >=0.0 && p<=1.0);
}

处理阶乘和大型乘法

任何大于 12! 的阶乘都会立即导致无符号整数溢出。所有概率类都涉及将一个非常大的数乘以一个非常小的数，以得到一个介于 0.0 和 1.0 之间的最终值。我利用这个知识在计算结果时将其优势发挥出来，从而避免溢出问题。类 ComputeResult() 方法中的每个公式都分解为其单独的因子。如果当前结果大于 1.0，那么我们知道需要除以某个数或乘以一个分数。如果当前结果小于 1.0，那么我们需要乘以一个因子。这种行为会导致结果在计算过程中在 1.0 值附近浮动。这也能减少有效数字的损失。大多数概率类都有一个需要计算的组合，它基本上是两个其他阶乘的商。一个有趣的且非常有用的数字抵消发生在组合计算中，这大大减少了乘法的次数。许多阶乘的使用通常需要除以其他阶乘或乘以小数，因此这种查看整个方程然后根据结果范围进行优化的通用方法是计算阶乘或其他通常具有非常大或非常小的中间结果的方程的好方法。

double BinomialProbability::ComputeResult() throw (ProbabilityException)
{
    if(m_trials == 0)
        return 0.0;
    // initialize some variables
    double result = 1.0;
    int Y = m_successes;
    int N = m_trials;
    double P = m_chance_of_success;
    double Q = 1.0 - P;
    int range = 0, np =0, nq = 0, nnumer = 0, ndenom = 0;
    // validate
    assert( Y<=N && Y >=0);
    assert( P <= 1.0 && P >=0.0);
    // check optimizations
    if(Y == 0){
        return result = pow(Q,N);
    } 
    if(Y == N){
        return result = pow(P,Y);
    }
    // reorder the factorials to account for cancellations
    // in numerator and denominator.
    if(Y < N-Y){
        range = Y;    // N-Y! cancels out
    }
    else{
        range = N-Y;    // Y! cancels out
    }
    np = Y;
    nq = N-Y;
    ndenom = range;
    nnumer = N;
    
    while(np > 0 || nq > 0 || ndenom > 0 || nnumer >(N-range)){
        // If the result is greater than one, we want to divide by 
        // a denominator digit or multiply by percentage p or q.
        // If we are out of numerator digits then finish multiplying
        // with our powers of p or q or dividing by a denom digit.
        if(result >= 1.0 || nnumer ==(N-range)){
            if(ndenom > 0){
                //m_resut *= (1.0/ndenom);
                result /= ndenom;
                --ndenom;
            }
            else if(nq > 0){
                result *= Q;
                --nq;
            }
            else if(np > 0){
                result *= P;
                --np;
            }
            else {
                throw(ProbabilityException(
                "Binomial Probability computation error- \
                check success percentage between 0 and 1"));
            } 
        }
        // If the result is less than one then we want to multiply
        // by a numerator digit. If we are out of denominator digits,
        // powers of p or powers of q then multiply rest of result 
        // by numerator digits.
        else if(result < 1.0 || np==0  ){
            if(nnumer >(N-range)){
                result *= nnumer;
                --nnumer;
            }
            else{
                throw(ProbabilityException(
                "Binomial Probability computation error- \ 
                unknown error"));
            }
        }
        else{
            throw(ProbabilityException(
            "Binomial Probability computation error- \
             possible value infinity or NaN"));
        }
    }
    return result; 
    
}

蒙特卡洛方法

最直观的概率计算方法之一是使用随机数生成器，然后计算命中或未命中的次数。再次尝试计算五次试验恰好掷出三个 1 的概率。考虑以下代码。

int Accuracy = 1000;
int N = 5;
int Y = 3;
int successess = 0;
for(int i =0; i < Accuracy; i++)
{
    int hits = 0;
    for(int j =0;j < N ; j++)
    {
        die = rand()%6 +1;  // give a random number from 1 to 6
        if(die == 1) //check for a one rolled
            hits++;
    }
    if(hits == Y) // check to see if 3 hits
        successes++;
}
double result = 1.0 * successes/Accuracy;

要获得此实验的百分比几率，您必须多次运行掷 5 个骰子的实验，然后计算掷出 3 个 1 的次数，再除以实验的总次数。此方法的问题在于没有一致性的保证。如果随机数生成器不像您期望的那样随机，那么准确性就无从谈起。要提高此方法的准确性，您必须多次运行实验。通常，为了获得不错的结果，您至少需要运行实验 1000 次。要获得更多小数位的准确性，您需要增加 Accuracy 的值。在复杂的情况下，所有这些都需要大量的 CPU 周期，精度越高，所需的周期就越多。

几何概率

此类概率模型用于计算事件成功（或失败）的试验次数。假设武器在比赛中有 2% 的失误率。那么在第 4 次射击时武器失误的概率是多少？这里我们只关心第 4 次射击，P(Y=4)，但我们也可以使用累积概率 P(Y<=4) 来查看在第 5 次射击之前发生失误的概率。当我们说 P(Y=4) 时，我们是在说“连续三次射击成功，然后一次失误的概率是多少？”对于我们的例子，参数很简单

    GeometricProbability gp(4, .02);
    // and for cumulative
    GeometricProbability gpc(4, .02, 
                Probability::RandomVariableType::LESSTHANOREQUAL);

负二项概率

此类概率与几何概率类似，模型计算的是第 K 个事件成功（或失败）的试验次数。就像上一个例子一样，我们可以问在第 20 次射击时武器失误三次的概率是多少。这意味着在前 19 次射击中有两次失误，而最后一次射击是第三次失误。

    NegativeBinomialProbability np(20,3,.02);

超几何概率

此类模型从两组中进行选择的更复杂的事件。它也是模拟现代彩票的概率。假设我们有一个袋子，里面有 3 个红球和 5 个黑球。我们想知道在抓取 3 次的情况下，选出 1 个红球的概率。袋子是黑的，并且假定抓取是随机的。超几何概率考虑 4 个参数。

N - 总体大小，共 8 个球，3 个红球和 5 个黑球
n - 样本大小，我们抓取 3 次（这假设抓取的球不会放回袋子！！）
r - 我们感兴趣的集合。我们想要红球，该集合中有 3 个。
y - 我们感兴趣的样本集中的项目数。我们只关心恰好得到一个。如果我们想考虑 1 个或更多，或 2 个或更少，或 3 个或更少，那么它将是累积概率。

   HyperGeometricProbability hgp(8,3,3,1);

一个更实际的例子可能是找出法律案件中某个事件的概率。从 20 名工人中，其中有三名是少数族裔，被选拔担任两种类型的职位，16 个普通劳动岗位和 4 个重体力劳动岗位。少数族裔在诉讼中声称所有人都被选为重体力劳动岗位，而雇主则认为分配是随机的。要计算所有少数族裔都被选入重体力劳动岗位的几率，我们将总体大小设为 N=20，样本大小设为 n=4（表示重体力劳动岗位的数量），我们感兴趣的集合是 r=3（表示少数族裔的数量），我们关心的是当 y=3 时，即所有少数族裔都被选中的情况。

   HyperGeometricProbability hgp(20,4,3,3);

泊松概率

泊松概率模型基于已知平均值来预测事件。平均值基于可以在任何时间发生且被认为是瞬时的事件（意味着两个事件不会在**同一时刻**发生）。例如，每页的打字错误数量或一年中街道交叉口的交通事故数量的平均值具有泊松“ish”的特性。如果我们知道某个交叉口平均每年发生 7 起交通事故，并且我们想输出今年发生 5 起或更少交通事故的概率……

cout<< PoissonProbability(5,7,
        Probability::RandomVariableType::LESSTHANOREQUAL).GetResult();

集合论与复合事件

这些概率类中未包含的一点是混合概率以获得自己相关答案的集合论。要模拟更复杂的事件，您需要理解子事件之间的独立性和互斥性概念。概率有关于事件之间独立性和互斥性的加法和乘法规则。这些涉及到集合论中的并集和交集概念。详细解释超出了本文的范围。读者还应该学习贝叶斯定理。网上有大量的示例和解释。
例如，用两个骰子掷出七的概率需要更多的思考。为了展示复合事件的模型，我们将需要将随机变量 Y1 和 Y2 的概率相乘，其中 Y1 + Y2 = 7，因为每次掷骰子的事件彼此独立。表示法为 P(Y1 = a) ∩ P(Y2 = 7-a)，即 P(Y1=a) * P(Y2=7-a)。

int Y1, Y2;
int N = 6;
int Roll = 7;
double p = 1.0/6.0;
double result = 0.0;
for(Y1=1; Y1<=6; Y1++)
{
    Y2 = Roll - Y1;
    BinomialProbability bp1( N, Y1, p);
    BinomialProbability bp2( N, Y2, p);
    // The answer is all the ways a seven can be 
    // rolled so we accumulate each independent
    // probability of P(y1=a)*P(y2=7-a)
    result += bp1.GetResult() * bp2.GetResult()
}
cout << result;

此外，每种概率都包含分布的期望值或均值以及方差。可以通过 GetExpectedValue 和 GetVariance 函数检索这些值。