本文链接：https://blog.csdn.net/hyq458941968/article/details/100055618

之前在看pjmedia的jbuf时，发现里面用到了pjlib中的数学统计，来实现抖动的测量。进去查看，发现方差的实现跟定义完全不一样，非常简洁，推导了1个多小时，终于推导出里面的代码实现原理。果然，程序员都是给数学家打工的。

pjlib实现的数学和统计文件是math.h，没有math.c文件，所有实现都内联在头文件。其中统计结构体如下：

/**
 * This structure describes statistics state.
 */
typedef struct pj_math_stat
{
    int		     n;		/* number of samples	*/    //数据集个数
    int		     max;	/* maximum value	*/        //数据集最大值
    int		     min;	/* minimum value	*/        //数据集最小值
    int		     last;	/* last value		*/        //最新的数据
    int		     mean;	/* mean			*/            //平均数整型

    /* Private members */
#if PJ_HAS_FLOATING_POINT
    float	     fmean_;	/* mean(floating point) */    //平均数浮点型
#else
    int		     mean_res_;	/* mean residu		*/
#endif
    pj_highprec_t    m2_;	/* variance * n		*/        //方差乘以个数
} pj_math_stat;

每次输入一个新数据，该数据集就能统计出个数，最大值最小值等，其中最重要的是平均值和方差乘以n的算法。但是代码中平均值和n倍方差不是按照定义计算的，而是经过化简。其中我们使用浮点型，PJ_HAS_FLOATING_POINT为1。

/**
 * Update statistics state as a new sample comes.
 *
 * @param stat	    Statistic state.
 * @param val	    The new sample data.
 */
PJ_INLINE(void) pj_math_stat_update(pj_math_stat *stat, int val)
{
#if PJ_HAS_FLOATING_POINT
    float	     delta;
#else
    int		     delta;
#endif

    stat->last = val;
    
    if (stat->n++) {
	if (stat->min > val)
	    stat->min = val;
	if (stat->max < val)
	    stat->max = val;
    } else {
	stat->min = stat->max = val;
    }

#if PJ_HAS_FLOATING_POINT
    delta = val - stat->fmean_;            //最新值减去上一次平均值
    stat->fmean_ += delta/stat->n;        //新的平均值=上一次平均值+(最新值-上一次平均值)/n
    
    /* Return mean value with 'rounding' */
    stat->mean = (int) (stat->fmean_ + 0.5);//整型值进1

    stat->m2_ += (int)(delta * (val-stat->fmean_));//最新方差值=上一次方差值+(最新值-上一次平均值)x(最新值-本次平均值)
#else
    delta = val - stat->mean;
    stat->mean += delta/stat->n;
    stat->mean_res_ += delta % stat->n;
    if (stat->mean_res_ >= stat->n) {
	++stat->mean;
	stat->mean_res_ -= stat->n;
    } else if (stat->mean_res_ <= -stat->n) {
	--stat->mean;
	stat->mean_res_ += stat->n;
    }

    stat->m2_ += delta * (val-stat->mean);
#endif
}

下面通过数学化解来解释上面的注释的代码

按定义，平均值的计算应该是 $\overline{x}=\frac{x_{1}+x_{2}+...+x_{n}}{n}=\frac{\sum_{1}^{n}x_{i}}{n}$ ，经过化简，平均值可采用差分方程的写法

从公式来看，新的平均值可以根据上一次的平均值求出来，而不必根据定义。

同样，方差的定义是 $s^{2}=\frac{ (x_{1}-\overline{x})^{2}+(x_{2}-\overline{x})^{2} +... +(x_{n}-\overline{x})^{2} }{n}=\frac{ \sum_{1}^{n}(x_{i}-\overline{x})^{2} }{n}$