logo标签列表

C语言中浮点数精度运算的问题

cfloating-pointmachine-learningspam-prevention
15

我正在为我的课程项目之一，用C语言实现“朴素贝叶斯分类器”。我的项目是使用大量的训练数据实现一个文档分类应用程序（特别是垃圾邮件）。

现在我遇到了问题，因为C语言的数据类型限制，无法实现该算法。

( 我使用的算法可以在这里找到：http://en.wikipedia.org/wiki/Bayesian_spam_filtering )

问题陈述： 该算法涉及对文档中的每个单词进行处理，并计算其作为垃圾单词的概率。如果P1、P2、P3....PN是单词1、2、3......N的概率。则根据下式计算出文档是否为垃圾邮件的概率：

alt text

在这里，概率值很容易接近0.01。所以即使我使用“double”数据类型，我的计算也会失败。为了确认这一点，我编写了下面的示例代码。

#define PROBABILITY_OF_UNLIKELY_SPAM_WORD     (0.01)
#define PROBABILITY_OF_MOSTLY_SPAM_WORD     (0.99)

int main()
{
    int index;
    long double numerator = 1.0;
    long double denom1 = 1.0, denom2 = 1.0;
    long double doc_spam_prob;

    /* Simulating FEW unlikely spam words  */
    for(index = 0; index < 162; index++)
    {
        numerator = numerator*(long double)PROBABILITY_OF_UNLIKELY_SPAM_WORD;
        denom2    = denom2*(long double)PROBABILITY_OF_UNLIKELY_SPAM_WORD;
        denom1    = denom1*(long double)(1 - PROBABILITY_OF_UNLIKELY_SPAM_WORD);
    }
    /* Simulating lot of mostly definite spam words  */
    for (index = 0; index < 1000; index++)
    {
        numerator = numerator*(long double)PROBABILITY_OF_MOSTLY_SPAM_WORD;
        denom2    = denom2*(long double)PROBABILITY_OF_MOSTLY_SPAM_WORD;
        denom1    = denom1*(long double)(1- PROBABILITY_OF_MOSTLY_SPAM_WORD);
    }
    doc_spam_prob= (numerator/(denom1+denom2));
    return 0;
}

我尝试了Float、double甚至是long double数据类型，但仍然存在同样的问题。因此，假设在我正在分析的10万个单词的文档中，只有162个单词具有1%的垃圾邮件概率，而其余99838个单词明显是垃圾邮件，则由于精度误差（因为分子很容易变为零），我的应用程序仍将其视为非垃圾邮件！

这是我第一次遇到这样的问题。那么这个问题应该如何解决呢？

6个回答
19

在机器学习中常常会出现这种情况。据我所知，无法避免精度损失。因此，我们使用log函数，并将除法和乘法转换为减法和加法来绕过这个问题。

于是我开始算了一下。

原始方程为：

Problem

我稍微修改了一下：

enter image description here

两边取对数：

enter image description here

设

enter image description here

代入得：

enter image description here

因此，计算联合概率的备选公式为：

enter image description here

如果需要我进一步解释，请留言。

3
在对概率进行计算时，使用对数域是唯一明智的方式。对数似然比（Jacob答案中的倒数第二个方程）是最容易处理的形式。 - Adam Bowen
2
@Microkernel：谢谢 :) - 你可以使用 exp 函数 http://www.codecogs.com/reference/c/math.h/exp.php 。即使用 exp(eta) 代替 pow(2.71828182845905,eta)。 - Jacob
6
当个体p_i很小时，这仍会失去一些准确性；如果这对你的目的有影响，解决方法是将ln(1-p_i)替换为log1p(-p_i)，这样就不会遇到相同的问题。（log1p是C标准库中未充分利用的宝石） - Stephen Canon
2
如果p_i的二进制指数是-n，在计算log（1-p_i）时，您应该期望失去n-1位精度。因此，如果p_i是0.1（二进制指数：-3），则您将失去2位精度，而如果使用log1p（-p_i）可能会保留更多精度。显然，这并不太糟糕，但如果p_i比这小得多，则损失可能相当大。是否值得担心这种差异取决于p_i的分布情况。如果它们都很小且规模相似，则问题非常重要。如果它们具有极大的规模差异，则可能完全没有关系。 - Stephen Canon
3
请注意，具体使用情况下这并不重要，因为当 p_i 很小时，log(p_i) 项将支配 log(1 - p_i) 项，所以小项的精度损失对最终结果影响微乎其微。在更一般的情况下，如果您有一个涉及形式为 log(1 + x) 的数值敏感计算，应考虑用 log1p(x) 替换它。 - Stephen Canon
显示剩余8条评论
4

这里有一个技巧：

for the sake of readability, let S := p_1 * ... * p_n and H := (1-p_1) * ... * (1-p_n), 
then we have:

  p = S / (S + H)
  p = 1 / ((S + H) / S)
  p = 1 / (1 + H / S)

let`s expand again:

  p = 1 / (1 +  ((1-p_1) * ... * (1-p_n)) / (p_1 * ... * p_n))
  p = 1 / (1 + (1-p_1)/p_1 * ... * (1-p_n)/p_n)

基本上，您将获得相当大的数字产品（介于 0 和对于 p_i = 0.01 时的 99 之间）。这个想法不是将大量小数互相乘以获得结果的 0，而是制作两个小数的商。例如，如果 n = 1000000 并且对于所有 i，p_i = 0.5，上述方法会给您 0/(0+0)，这是 NaN，而建议的方法会给您 1/(1+1*...1)，即 0.5。
当所有 p_i 都排序并且您按相反顺序配对它们时（假设 p_1 < ...<p_n），则以下公式将获得更好的精度：您可以获得更好的结果。
  p = 1 / (1 + (1-p_1)/p_n * ... * (1-p_n)/p_1)

这样你就可以将大分子（小的p_i）与大分母（大的p_(n+1-i)）相除，将小分子与小分母相除。

编辑：MSalter在他的回答中提出了一个有用的进一步优化。使用它，公式如下：

  p = 1 / (1 + (1-p_1)/p_n * (1-p_2)/p_(n-1) * ... * (1-p_(n-1))/p_2 * (1-p_n)/p_1)
这是一个非常有趣的想法... 我会尝试一下并查看Jacob的回答，以确定哪一个更符合我的要求。非常感谢 :) - Microkernel
“排序术”确实有效，但如果您动态选择大或小的术语来保持中间结果约为1.0，则效果更好。请参阅我的答案。 - MSalters
@MSalters：说得好。我认为最好的解决方案是像我做的那样，将概率成对相反地排列，以保持因子更接近1，然后按照你提出的方式交替重新排列因子。 - back2dos
实际上，我最初也采用了同样的方法，但后来发现如果你有一小部分极端项和大量非极端项，就会出现失控效应。即少量p=0.01与大量p=0.51相平衡。最初，你会将少量的0.01项与0.51项配对，并向无限远处奔去。之后，你会将那些p=0.51项配对，并不断地将无限乘以0.98。这根本行不通。 - MSalters
3
你的问题是因为你收集了太多没有考虑大小的术语。一种解决方法是取对数，另一种方法是对单个术语进行排序。首先，让我们将方程重写为1/p = 1 + ∏((1-p_i)/p_i)。现在你的问题是有些术语很小，而其他一些术语很大。如果你有太多相同大小的术语连续在一起，就会下溢，在有太多大的术语时你会溢出中间结果。
所以，不要把太多相同大小的术语放在一起。将术语(1-p_i)/p_i进行排序。结果，第一个术语是最小的术语，最后一个是最大的术语。现在，如果你立即将它们相乘，仍然会下溢。但计算顺序并不重要。使用两个迭代器进入临时集合。一个从开始位置开始（即(1-p_0)/p_0），另一个从结束位置开始（即(1-p_n)/p_n），你的中间结果从1.0开始。现在，当你的中间结果>=1.0时，从前面取一个术语，当你的中间结果<1.0时，从后面取一个结果。
结果是，随着你取的术语，中间结果将在1.0周围振荡。当你用完小的或大的术语时，它只会上升或下降。但那没关系。此时，你已经消耗了两端的极值，所以中间结果会慢慢接近最终结果。
当然也有真正的溢出可能性。如果输入完全不太可能是垃圾邮件（p=1E-1000），那么1/p会溢出，因为∏((1-p_i)/p_i)会溢出。但由于术语已排序，我们知道中间结果仅在∏((1-p_i)/p_i)溢出时会溢出。所以，如果中间结果溢出，就没有后续的精度损失。
+1. 我更新了我的答案。我认为最好的方法是将两种算法结合起来，因为我的算法在计算因子时精度损失较小，而你的算法在计算整体乘积时精度损失较小。 - back2dos
2

尝试计算逆数1/p。这将给您一个形式为1 + 1/(1-p1)*(1-p2)...的方程。

如果您计算每个概率的出现次数--看起来您有一小部分重复值--您可以使用pow()函数--pow(1-p, p的出现次数)*pow(1-q, q的出现次数)--并避免每次乘法的单独舍入误差。

那不是1/p，看看我的答案。即使你是对的，它仍然涉及乘以(1-p_i)，它可以取任何从0到1的值，所以如果它取值接近1，我们又回到了起点。 - Jacob
1

您可以使用百分数或千分数的概率：

doc_spam_prob= (numerator*100/(denom1+denom2));

或者

doc_spam_prob= (numerator*1000/(denom1+denom2));

或者使用其他系数

0
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接