使用Boost在C++中计算样本向量的平均值和标准差

Question

有没有一种方法可以使用Boost计算包含样本的向量的平均值和标准差？

还是我必须创建一个累加器并将向量馈入其中？

- user393144

Answer 1

我不知道Boost库是否有更具体的函数，但您可以使用标准库来实现。

给定std::vector<double> v，以下是一种简单的方法：

#include <numeric>

double sum = std::accumulate(v.begin(), v.end(), 0.0);
double mean = sum / v.size();

double sq_sum = std::inner_product(v.begin(), v.end(), v.begin(), 0.0);
double stdev = std::sqrt(sq_sum / v.size() - mean * mean);

当值过大或过小时，这种方法容易发生溢出或下溢。稍微更好的计算标准差的方法是：

double sum = std::accumulate(v.begin(), v.end(), 0.0);
double mean = sum / v.size();

std::vector<double> diff(v.size());
std::transform(v.begin(), v.end(), diff.begin(),
               std::bind2nd(std::minus<double>(), mean));
double sq_sum = std::inner_product(diff.begin(), diff.end(), diff.begin(), 0.0);
double stdev = std::sqrt(sq_sum / v.size());

C++11的更新内容:

调用std::transform可以使用lambda函数来替代std::minus和std::bind2nd(现已弃用):

std::transform(v.begin(), v.end(), diff.begin(), [mean](double x) { return x - mean; });

- musiphil

---

第一组方程式不起作用。我输入了10和2，得到的输出是4。乍一看，我认为这是因为它假定(a-b)^2 = a^2-b^2。 - Charles L.

@CharlesL.：它应该可以工作，4是正确答案。 - musiphil

@StudentT：不行，但你可以在上面的最后一行中将 v.size() 替换为 (v.size() - 1)：std::sqrt(sq_sum / (v.size() - 1))。（对于第一种方法，有一点复杂：std::sqrt(sq_sum / (v.size() - 1) - mean * mean * v.size() / (v.size() - 1))。） - musiphil

使用std::inner_product计算平方和非常巧妙。 - Paul R

我可以第一手证实，第一个实现对于极小的数字确实会发生溢出/下溢。我不得不改用第二个实现，然后标准差就没有得到NAN值了。为了避免溢出/下溢，这两行额外的代码是值得的！ - Trevor Boyd Smith

Answer 2

如果性能对你很重要，并且你的编译器支持lambda，stdev计算可以更快、更简单地完成：在使用VS 2012进行测试时，我发现以下代码比所选答案中给出的Boost代码快10倍以上；它也比musiphil提供的使用标准库的更安全的答案快5倍以上。

注意，我正在使用样本标准差，因此下面的代码给出稍微不同的结果（为什么标准偏差中有一个负号）。

double sum = std::accumulate(std::begin(v), std::end(v), 0.0);
double m =  sum / v.size();

double accum = 0.0;
std::for_each (std::begin(v), std::end(v), [&](const double d) {
    accum += (d - m) * (d - m);
});

double stdev = sqrt(accum / (v.size()-1));

- Josh Greifer

---

感谢您即使一年后仍然分享这个答案。现在我又来了一年，将其变成了通用的值类型和容器类型。请看这里（注意：我猜我的基于范围的for循环与您的lambda代码一样快。） - leemes

使用std::end(v)和v.end()有什么区别？ - spurra

std::end()函数是C++11标准新增的，用于处理没有v.end()这样的情况。对于较不常见的容器，std::end可以进行重载--请参见http://en.cppreference.com/w/cpp/iterator/end。 - pepr

我实际上是在尝试理解为什么它比“更安全”的答案快5倍。你的答案易于理解，而且凭借我有限的统计知识，我知道它会起作用。我知道时间复杂度是相同的。当你说比musiphil的更安全的答案更快时，你是指他的非朴素（第二个）答案吗？还是比朴素（第一个）答案更快？ - dev_nut

首先，“安全”的答案（就像我的答案一样）需要对数组进行三次遍历：一次求和，一次求差的平均值，和一次平方。而在我的代码中，只需要两次遍历--将这两次遍历合并为一次。而且（截至我上次查看时，已经有相当长一段时间了！），inner_product调用没有被优化掉。此外，“安全”代码将v复制到一个全新的diffs数组中，这会增加更多的延迟。我认为我的代码更易读，而且容易移植到JavaScript和其他语言中 :) - Josh Greifer

Answer 3

在Boost中，使用累加器是计算平均值和标准差的方法。

accumulator_set<double, stats<tag::variance> > acc;
for_each(a_vec.begin(), a_vec.end(), bind<void>(ref(acc), _1));

cout << mean(acc) << endl;
cout << sqrt(variance(acc)) << endl;

 

- David Nehme

---

注意，tag::variance使用近似公式计算方差。tag::variance(lazy)使用精确公式计算，即“二阶矩减平均数的平方”，如果由于舍入误差导致方差非常小，则会产生不正确的结果。它实际上可以产生负方差。 - panda-34

如果你知道你将要处理大量数字，那么请使用递归（在线）算法。这将解决下溢和上溢问题。 - Kemin Zhou

Answer 4

优化 musiphil的答案后，您可以编写一个标准差函数，而无需使用临时向量diff。只需使用C++11的lambda功能，通过单个inner_product调用即可实现：

double stddev(std::vector<double> const & func)
{
    double mean = std::accumulate(func.begin(), func.end(), 0.0) / func.size();
    double sq_sum = std::inner_product(func.begin(), func.end(), func.begin(), 0.0,
        [](double const & x, double const & y) { return x + y; },
        [mean](double const & x, double const & y) { return (x - mean)*(y - mean); });
    return std::sqrt(sq_sum / func.size() - 1);
}

我怀疑多次执行减法比使用额外的中间存储更便宜，而且我认为这样更易读，但我还没有测试过性能。

至于为什么要使用 N-1（如 func.size() - 1），请参阅这些 问题 - 注意问题中提到我们有一个“包含样本”的向量。

- codeling

---

我认为这个计算得到的是方差，而不是标准差。 - sg_man

标准差是通过除以N而不是N-1来计算的。 为什么要将sq_sum除以func.size()-1？ - pocjoc

为什么是N-1？在计算总体标准差时，你要除以N，而在计算样本标准差时，你要除以N-1。https://stats.stackexchange.com/q/3931 https://stats.stackexchange.com/q/485326 - Jonathon Reinhart

@JonathonReinhart 你的意思是说使用 N-1 进行初始计算是正确的，对吗？ - codeling

我在说，这两者都是正确的，取决于你正在计算什么。 - Jonathon Reinhart

Answer 5

看起来下面这个优美的递归解决方案还没有被提及，尽管它已经存在了很长时间。参考Knuth的《计算机程序设计艺术》。

mean_1 = x_1, variance_1 = 0;            //initial conditions; edge case;

//for k >= 2, 
mean_k     = mean_k-1 + (x_k - mean_k-1) / k;
variance_k = variance_k-1 + (x_k - mean_k-1) * (x_k - mean_k);

对于一个列表中含有 n≥2 个值的情况，标准差的估计值为：

stddev = std::sqrt(variance_n / (n-1)). 

希望这可以帮到你！

- galactica

---

这很酷。我用索引循环实现了它（https://pastebin.com/aRd1ChjD），但它的运行速度比基于stl的解决方案慢三倍。 - DarioP

Answer 6

我和Josh Greifer的答案类似，但是更为普遍，适用于样本协方差。 样本方差只是样本协方差的两个输入相同而已。这包括Bessel的相关性。

    template <class Iter> typename Iter::value_type cov(const Iter &x, const Iter &y)
    {
        double sum_x = std::accumulate(std::begin(x), std::end(x), 0.0);
        double sum_y = std::accumulate(std::begin(y), std::end(y), 0.0);

        double mx =  sum_x / x.size();
        double my =  sum_y / y.size();

        double accum = 0.0;

        for (auto i = 0; i < x.size(); i++)
        {
            accum += (x.at(i) - mx) * (y.at(i) - my);
        }

        return accum / (x.size() - 1);
    }

- ABCD

Answer 7

比之前提到的版本快2倍 - 主要是因为将transform()和inner_product()循环合并了。对于我的快捷方式/typedefs/macro感到抱歉：Flo = 浮点数，CR const ref，VFlo-向量。已在VS2010中进行了测试。

#define fe(EL, CONTAINER)   for each (auto EL in CONTAINER)  //VS2010
Flo stdDev(VFlo CR crVec) {
    SZ  n = crVec.size();               if (n < 2) return 0.0f;
    Flo fSqSum = 0.0f, fSum = 0.0f;
    fe(f, crVec) fSqSum += f * f;       // EDIT: was Cit(VFlo, crVec) {
    fe(f, crVec) fSum   += f;
    Flo fSumSq      = fSum * fSum;
    Flo fSumSqDivN  = fSumSq / n;
    Flo fSubSqSum   = fSqSum - fSumSqDivN;
    Flo fPreSqrt    = fSubSqSum / (n - 1);
    return sqrt(fPreSqrt);
}

- slyy2048

---

Cit()循环可以写成以下形式吗？for(float f: crVec) { fSqSum += f * f; fSum += f; } - Elfen Dew

在C++11中是可以的。 尝试使用使其版本无关的宏。 更新了代码。 PS. 为了可读性，我通常更喜欢每个LOC执行1个操作。编译器应该看到这些是常量迭代，并在它“认为”一次迭代更快时将它们合并。 通过小而短的步骤来完成（例如不使用std :: inner_product（）），类似于汇编风格，向新读者解释其含义。二进制文件将因此而变小（在某些情况下）。 - slyy2048

尝试使用使其版本无关的宏，但您仍限制自己使用非标准的Visual C++“for each”结构（https://dev59.com/DnVC5IYBdhLWcg3wvT5a） - codeling

@codeling 这只是一个宏，用于说明该帖子中的C++版本。那是算法 - 不是编码标准。当时我甚至使用了更丑陋的Cit(CFlo, crVec)，其中默认const-iter“cit”，但重新指示容器类型。 当可移植性成问题时，列出所有编译器/操作系统特定的宏是很好的选择。在使用boost的示例中，将其移植到std C ++也不容易。 我没有解释过短的Flo、VFlo、CR、SZ -> float、vector<float>、const&、size - 为了缩短std C++的迭代行。 同样的风格Crit(MSZPFlo, crMap) foo(*crit.second); //rev-iter - slyy2048

Answer 8

为了更精确地计算样本均值，可以使用以下r步递归：

mean_k=1/k*[(k-r)*mean_(k-r) + sum_over_i_from_(n-r+1)_to_n(x_i)],

其中选择r以使求和分量彼此更接近。

- Sergey Fedortsov

Answer 9

创建自己的容器：

template <class T>
class statList : public std::list<T>
{
    public:
        statList() : std::list<T>::list() {}
        ~statList() {}
        T mean() {
           return accumulate(begin(),end(),0.0)/size();
        }
        T stddev() {
           T diff_sum = 0;
           T m = mean();
           for(iterator it= begin(); it != end(); ++it)
               diff_sum += ((*it - m)*(*it -m));
           return diff_sum/size();
        }
};

它确实有一些限制，但是当您知道自己在做什么时，它的表现非常出色。

- Sushant Kondguli

---

回答这个问题的原因是：因为完全没有必要。与编写自由函数相比，创建自己的容器绝对没有任何好处。 - Konrad Rudolph

我甚至不知道该从哪里开始。您正在使用列表作为底层数据结构，甚至没有缓存值，这是我能想到的使用类似容器结构的少数原因之一。特别是如果值很少改变且平均值/标准差经常需要。 - Creat

Answer 10

//在C++中表示偏差

//观测值与感兴趣数量（如总体均值）的真实值之间的差异是误差，观测值与真实值估计值（例如样本均值）之间的差异是残差。这些概念适用于区间和比例测量水平的数据。

#include <iostream>
#include <conio.h>
using namespace std;

/* run this program using the console pauser or add your own getch,     system("pause") or input loop */

int main(int argc, char** argv)
{
int i,cnt;
cout<<"please inter count:\t";
cin>>cnt;
float *num=new float [cnt];
float   *s=new float [cnt];
float sum=0,ave,M,M_D;

for(i=0;i<cnt;i++)
{
    cin>>num[i];
    sum+=num[i];    
}
ave=sum/cnt;
for(i=0;i<cnt;i++)
{
s[i]=ave-num[i];    
if(s[i]<0)
{
s[i]=s[i]*(-1); 
}
cout<<"\n|ave - number| = "<<s[i];  
M+=s[i];    
}
M_D=M/cnt;
cout<<"\n\n Average:             "<<ave;
cout<<"\n M.D(Mean Deviation): "<<M_D;
getch();
return 0;

}

- ali