练习43：一个简单的统计引擎

原文：Exercise 43: A Simple Statistics Engine

这是一个简单的算法，我将其用于“联机”（不储存任何样本）收集概要统计。我在任何需要执行一些统计，比如均值、标准差和求和中使用它，但是其中我并不会储存所需的全部样本。我只需要储存计算出的结果，它们仅仅含有5个数值。

计算标准差和均值

首先你需要一系列样本。它可以使任何事情，比如完成一个任务所需的时间，某人访问某个东西的次数，或者甚至是网站的评分。是什么并不重要，只要你能得到一些数字，并且你想要知道它们的下列概要统计值：

sum

对所有数字求和。

sumsq（平方和）

对所有数字求平方和。

count(n)

求出样本数量。

min

求出样本最小值。

max

求出样本最大值。

mean

求出样本的均值。它类似于但又不是中位数，但可作为中位数的估计。

stddev

使用$sqrt(sumsq - (sum * mean) / (n - 1) )来计算标准差，其中sqrt为math.h头文件中的平方根。

我将会使用R来验证这些计算，因为我知道R能够计算正确。

> s <- runif(n=10, max=10)
> s
 [1] 6.1061334 9.6783204 1.2747090 8.2395131 0.3333483 6.9755066 1.0626275
 [8] 7.6587523 4.9382973 9.5788115
> summary(s)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
 0.3333  2.1910  6.5410  5.5850  8.0940  9.6780
> sd(s)
[1] 3.547868
> sum(s)
[1] 55.84602
> sum(s * s)
[1] 425.1641
> sum(s) * mean(s)
[1] 311.8778
> sum(s * s) - sum(s) * mean(s)
[1] 113.2863
> (sum(s * s) - sum(s) * mean(s)) / (length(s) - 1)
[1] 12.58737
> sqrt((sum(s * s) - sum(s) * mean(s)) / (length(s) - 1))
[1] 3.547868
>

你并不需要懂得R，只需要看着我拆分代码来解释如何检查这些运算：

lines 1-4

我使用runit函数来获得“随机形式”的数字分布，之后将它们打印出来。我会在接下来的单元测试中用到它。

lines 5-7

这个就是概要，便于你看到R如何计算它们。

lines 8-9

这是使用sd函数计算的stddev。

lines 10-11

现在我开始手动进行这一计算，首先计算sum。

lines 12-13

stddev公式中的下一部分是sumsq，我可以通过sum(s * s)来得到，它告诉R将整个s列表乘以其自身，之后计算它们的sum。R的可以在整个数据结构上做运算，就像这样。

lines 14-15

观察那个公式，我之后需要sum乘上mean，所以我执行了sum(s) * mean(s)。

lines 16-17

我接着将sumsq参与运算，得到sum(s * s) - sum(s) * mean(s)。

lines 18-19

还需要除以n - 1，所以我执行了(sum(s * s) - sum(s) * mean(s)) / (length(s) - 1)。

lines 20-21

随后，我使用sqrt算出平方根，并得到3.547868，它符合R通过sd的运算结果。

实现

这就是计算stddev的方法，现在我可以编写一些简单的代码来实现这一计算。

#ifndef lcthw_stats_h
#define lctwh_stats_h

typedef struct Stats {
    double sum;
    double sumsq;
    unsigned long n;
    double min;
    double max;
} Stats;

Stats *Stats_recreate(double sum, double sumsq, unsigned long n, double min, double max);

Stats *Stats_create();

double Stats_mean(Stats *st);

double Stats_stddev(Stats *st);

void Stats_sample(Stats *st, double s);

void Stats_dump(Stats *st);

#endif

这里你可以看到我将所需的统计量放入一个struct，并且创建了用于处理样本和获得数值的函数。实现它只是转换数字的一个练习：

#include <math.h>
#include <lcthw/stats.h>
#include <stdlib.h>
#include <lcthw/dbg.h>

Stats *Stats_recreate(double sum, double sumsq, unsigned long n, double min, double max)
{
    Stats *st = malloc(sizeof(Stats));
    check_mem(st);

    st->sum = sum;
    st->sumsq = sumsq;
    st->n = n;
    st->min = min;
    st->max = max;

    return st;

error:
    return NULL;
}

Stats *Stats_create()
{
    return Stats_recreate(0.0, 0.0, 0L, 0.0, 0.0);
}

double Stats_mean(Stats *st)
{
    return st->sum / st->n;
}

double Stats_stddev(Stats *st)
{
   return sqrt( (st->sumsq - ( st->sum * st->sum / st->n)) / (st->n - 1) );
}

void Stats_sample(Stats *st, double s)
{
    st->sum += s;
    st->sumsq += s * s;

    if(st->n == 0) {
        st->min = s;
        st->max = s;
    } else {
        if(st->min > s) st->min = s;
        if(st->max < s) st->max = s;
    }

    st->n += 1;
}

void Stats_dump(Stats *st)
{
    fprintf(stderr, "sum: %f, sumsq: %f, n: %ld, min: %f, max: %f, mean: %f, stddev: %f",
            st->sum, st->sumsq, st->n, st->min, st->max,
            Stats_mean(st), Stats_stddev(st));
}

下面是stats.c中每个函数的作用：

Stats_recreate

我希望从一些数据中加载这些数据，这和函数让我重新创建Stats结构体。

Stats_create

只是以全0的值调用Stats_recreate。

Stats_mean

使用sum和n计算均值。

Stats_stddev

实现我之前的公式，唯一的不同就是我使用t->sum / st->n来计算均值，而不是调用Stats_mean。

Stats_sample

它用于在Stats结构体中储存数值。当你向它提供数值时，它看到n是0，并且相应地设置min和max。之后的每次调用都会使sum、sumsq和n增加，并且计算出这一新的样本的min和max值。

Stats_dump

简单的调试函数，用于转储统计量，便于你看到它们。

我需要干的最后一件事，就是确保这些运算正确。我打算使用我的样本，以及来自于R会话中的计算结果创建单元测试，来确保我会得到正确的结果。

#include "minunit.h"
#include <lcthw/stats.h>
#include <math.h>

const int NUM_SAMPLES = 10;
double samples[] = {
    6.1061334, 9.6783204, 1.2747090, 8.2395131, 0.3333483,
    6.9755066, 1.0626275, 7.6587523, 4.9382973, 9.5788115
};

Stats expect = {
    .sumsq = 425.1641,
    .sum = 55.84602,
    .min = 0.333,
    .max = 9.678,
    .n = 10,
};
double expect_mean = 5.584602;
double expect_stddev = 3.547868;

#define EQ(X,Y,N) (round((X) * pow(10, N)) == round((Y) * pow(10, N)))

char *test_operations()
{
    int i = 0;
    Stats *st = Stats_create();
    mu_assert(st != NULL, "Failed to create stats.");

    for(i = 0; i < NUM_SAMPLES; i++) {
        Stats_sample(st, samples[i]);
    }

    Stats_dump(st);

    mu_assert(EQ(st->sumsq, expect.sumsq, 3), "sumsq not valid");
    mu_assert(EQ(st->sum, expect.sum, 3), "sum not valid");
    mu_assert(EQ(st->min, expect.min, 3), "min not valid");
    mu_assert(EQ(st->max, expect.max, 3), "max not valid");
    mu_assert(EQ(st->n, expect.n, 3), "max not valid");
    mu_assert(EQ(expect_mean, Stats_mean(st), 3), "mean not valid");
    mu_assert(EQ(expect_stddev, Stats_stddev(st), 3), "stddev not valid");

    return NULL;
}

char *test_recreate()
{
    Stats *st = Stats_recreate(expect.sum, expect.sumsq, expect.n, expect.min, expect.max);

    mu_assert(st->sum == expect.sum, "sum not equal");
    mu_assert(st->sumsq == expect.sumsq, "sumsq not equal");
    mu_assert(st->n == expect.n, "n not equal");
    mu_assert(st->min == expect.min, "min not equal");
    mu_assert(st->max == expect.max, "max not equal");
    mu_assert(EQ(expect_mean, Stats_mean(st), 3), "mean not valid");
    mu_assert(EQ(expect_stddev, Stats_stddev(st), 3), "stddev not valid");

    return NULL;
}

char *all_tests()
{
    mu_suite_start();

    mu_run_test(test_operations);
    mu_run_test(test_recreate);

    return NULL;
}

RUN_TESTS(all_tests);

这个单元测试中没什么新东西，除了EQ宏。我比较懒，并且不想查询比较两个double值的标准方法，所以我使用了这个宏。double的问题是等性不是完全相等，因为我使用了两个不同的系统，并带有不同的四舍五入的位数。解决方案就是判断两个数“乘以10的X次方是否相等”。

我使用EQ来计算数字的10的幂，之后使用round函数来获得证书。这是个简单的方法来四舍五入N位小数，并以整数比较结果。我确定有数以亿计的其它方法能做相同的事情，但是现在我就用这种。

预期结果储存在Stats struct中，之后我只是确保我得到的数值接近R给我的数值。

如何使用

你可以使用标准差和均值来决定一个新的样本是否是“有趣”的，或者你可以使用它们计算统计量的统计量。前者对于人们来说更容易理解，所以我用登录的例子来做个简短的解释。

假设你在跟踪人们花费多长时间在一台服务器上，并且你打算用统计来分析它。每次有人登录进来，你都对它们在这里的时长保持跟踪，之后调用Stats_sample函数。我会寻找停留“过长”时间的人，以及“过短”的人。

比起设定特殊的级别，我更倾向于将一个人的停留时间与mean (plus or minus) 2 * stddev这个范围进行比较。我计算出mean和2 * stddev，并且如果它们在这个范围之外，我就认为是“有趣”的。由于我使用了联机算法来维护这些统计量，所以它非常快，并且我可以使软件标记在这个范围外的用户。

这不仅仅用于找出行为异常的用户，更有助于标记一些潜在的问题，你可以查看它们来观察发生了什么。它基于所有用户的行为来计算，这也避免了你任意挑出一个数值而并不基于实际情况的问题。

你可以从中学到的通用规则是，mean (plus or minus) 2 * stddev是90%的值预期所属的范围预测值，任何在它之外的值都是有趣的。

第二种利用这些统计量的方式就是继续将其用于其它的Stats计算。基本上像通常一样使用Stats_sample，但是之后在min、max、n、mean和stddev上执行Stats_sample。这会提供二级的度量，并且让你对比样本的样本。

被搞晕了吗？我会以上面的例子基础，并且假设你拥有100台服务器，每台都运行一个应用。你已经在每个应用服务器上跟踪了用户的登录时长，但是你想要比较所有的这100和应用，并且标记它们当中任何登录时间过长的用户。最简单的方式就是每次有人登录进来时，计算新的登录统计量，之后将Stats structs的元素添加到第二个Stats中。

你最后应该会得到一些统计量，它们可以这样命名：

均值的均值

这是一个Stats struct，它向你提供所有服务器的均值的mean和stddev。你可以用全局视角来观察任何在此之外的用户或服务器。

标准差的均值

另一个Stats struct，计算这些服务器的分布的统计量。你之后可以分析每个服务器并且观察是否它们中的任何服务器具有异常分散的分布，通过将它们的stddev和这个mean of stddevs统计量进行对比。

你可以计算出全部统计量，但是这两个是最有用的。如果你打算监视服务器上的移除登录时间，你可以这样做：

• 用户John登录并登出服务器A。获取服务器A的统计量，并更新它们。
• 获取mean of means统计量，计算出A的均值并且将其加入样本。我叫它m_of_m。
• 获取mean of stddev统计量，将A的标准差添加到样本中。我叫它m_of_s。
• 如果A的mean在m_of_m.mean + 2 * m_of_m.stddev范围外，标记它可能存在问题。
• 如果A的stddev在m_of_s.mean + 2 * m_of_s.stddev范围外，标记它可能存在行为异常。
• 最后，如果John的登录时长在A的范围之外，或A的m_of_m范围之外，标记为有趣的。

通过计算“均值的均值”，或者“标准差的均值”，你可以以最小的执行和储存总量，有效地跟踪许多度量。

附加题

• 将Stats_stddev 和 Stats_mean转换为static inline函数，放到stats.h文件中，而不是stats.c文件。
• 使用这份代码来编写string_algos_test.c的性能测试。使它为可选的，并且运行基准测试作为一系列样本，之后报告结果。
• 编写它的另一个语言的版本。确保这个版本基于我的数据正确执行。
• 编写一个小型程序，它能从文件读取所有数字，并执行这些统计。
• 使程序接收一个数据表，其中第一行是表头，剩下的行含有任意数量空格分隔的数值。你的程序应该按照表头中的名称，打印出每一列的统计值。