关于原子哥的一段求方差代码的疑惑和讨论

jiang039

大家好，我最近在努力学习原子哥的MiniFly，在参考原子哥的一段代码的时候，遇到了些问题，希望大神们能够帮我解开迷津！
下面复制一下原子哥的代码：
//计算方差和平均值
static void sensorsCalculateVarianceAndMean(BiasObj*bias,Axis3f*varOut,Axis3f*meanOut)
{
        u32 i;
        int64_t        sum[3]={0};
        int64_t        sumsq[3]={0};
       
        for (i=0;i<SENSORS_NBR_OF_BIAS_SAMPLES;i++)
        {
                sum[0]+=bias->buffer.x;
                sum[01]+=bias->buffer.y;
                sum[02]+=bias->buffer.z;
                sumsq[0]+=bias->buffer.x*bias->buffer.x;
                sumsq[01]+=bias->buffer.y*bias->buffer.y;
                sumsq[02]+=bias->buffer.z*bias->buffer.z;
        }
       
        varOut->x=(sumsq[0]-((int64_t)sum[0]*sum[0])/SENSORS_NBR_OF_BIAS_SAMPLES);
        varOut->y=(sumsq[01]-((int64_t)sum[01]*sum[01])/SENSORS_NBR_OF_BIAS_SAMPLES);
        varOut->z=(sumsq[02]-((int64_t)sum[02]*sum[02])/SENSORS_NBR_OF_BIAS_SAMPLES);

        meanOut->x=(float)sum[0]/SENSORS_NBR_OF_BIAS_SAMPLES;
        meanOut->y=(float)sum[1]/SENSORS_NBR_OF_BIAS_SAMPLES;
        meanOut->z=(float)sum[2]/SENSORS_NBR_OF_BIAS_SAMPLES;
}


这其中的“varOut->x”应该就是方差的意思了，那么问题就来了，按照原子哥的代码我们可以理解成：方差=平方和-和平方/数据个数,但是经过我的反复推导方差公式我发现
：方差=（平方和-和平方/数据个数）*（1/数据个数），这不禁让我反思起来我们开始的判断“varOut->x”指的是方差吗？

xianshasaman

jiang039 发表于 2018-9-6 14:57
为什么我推导出来的公式需要多乘以一个"1/n"

不会多一个1/n，应该多一个1/（n-1）
因为n-1在程序里是个立即寻址的常数，而在局部也只是个独立的系数，所以可以在函数外面约掉。这个和直接用整数表示保留一位小数的浮点数是一个道理。

ssssssssssss

[mw_shl_code=c,true]static void sensorsCalculateVarianceAndMean(BiasObj*bias,Axis3f*varOut,Axis3f*meanOut)
{
        u32 i;
        int64_t        sum[3]={0};
        int64_t        sumsq[3]={0};
        
        for (i=0;i<SENSORS_NBR_OF_BIAS_SAMPLES;i++)
        {
                sum[0]+=bias->buffer.x;
                sum[01]+=bias->buffer.y;
                sum[02]+=bias->buffer.z;
                sumsq[0]+=bias->buffer.x*bias->buffer.x;
                sumsq[01]+=bias->buffer.y*bias->buffer.y;
                sumsq[02]+=bias->buffer.z*bias->buffer.z;
        }
        
        varOut->x=(sumsq[0]-((int64_t)sum[0]*sum[0])/SENSORS_NBR_OF_BIAS_SAMPLES);
        varOut->y=(sumsq[01]-((int64_t)sum[01]*sum[01])/SENSORS_NBR_OF_BIAS_SAMPLES);
        varOut->z=(sumsq[02]-((int64_t)sum[02]*sum[02])/SENSORS_NBR_OF_BIAS_SAMPLES);

        meanOut->x=(float)sum[0]/SENSORS_NBR_OF_BIAS_SAMPLES;
        meanOut->y=(float)sum[1]/SENSORS_NBR_OF_BIAS_SAMPLES;
        meanOut->z=(float)sum[2]/SENSORS_NBR_OF_BIAS_SAMPLES;
}[/mw_shl_code]

edmund1234

基于证据挑战权威， 正是科学精神之真缔， 赞一个
但我觉得原来的没错啊， 你的证明有点看不懂

xianshasaman

你看错了吧，c语言除法优先级比减法高，人家的代码是先算的除法。

还有，你对样本方差的理解也不太准确，样本方差公式是这样的：

注意外面的分母是n-1，变形后里面的分母是n，代码里面的 SENSORS_NBR_OF_BIAS_SAMPLES 这个宏就是n，n-1在函数外面被约掉了。
至于为什么是n-1，你可以参考下概率统计的教材里面无偏估计量的内容。

jiang039

edmund1234 发表于 2018-9-6 12:22
基于证据挑战权威， 正是科学精神之真缔， 赞一个
但我觉得原来的没错啊， 你的证明有点看不懂

谢谢大哥在精神上对我的赞赏，我理解的原子哥的公式是：用数据的平方和-（数据和的平方/数据个数）得到的方差，但是经过我反复推导，始终推导不出这个公式，

xianshasaman

补一个推到过程

jiang039

xianshasaman 发表于 2018-9-6 14:23
你看错了吧，c语言除法优先级比减法高，人家的代码是先算的除法。

还有，你对样本方差的理解也不太准确 ...

谢谢大哥的指点，我很仔细的看了您的解释，我知道c语言里是先算的除法，就拿这段代码来说，它是先算的除法后算的减法，简单的是S^2=sumsq-(sum^2/n),其中“sumsq”是平方和、“sum”是和。而我却推导出这样一个公式：S^2=（sumsq-(sum^2/n)）*1/n；所以我一直不是很理解

jiang039

jiang039 发表于 2018-9-6 14:52
谢谢大哥的指点，我很仔细的看了您的解释，我知道c语言里是先算的除法，就拿这段代码来说，它是先算的除 ...

为什么我推导出来的公式需要多乘以一个"1/n"

jiang039

xianshasaman 发表于 2018-9-6 14:50
补一个推到过程

大哥你这个我可以理解为：【样本的平方和-样本和的平方/样本个数】/自由度 对吧？
但是原子哥的代码里没有除掉自由度，也就是没有你推导的（1/n-1）,我推导的（1/n）,所以最终得出来的结果还是方差吗？

jiang039

xianshasaman 发表于 2018-9-6 15:05
不会多一个1/n，应该多一个1/（n-1）
因为n-1在程序里是个立即寻址的常数，而在局部也只是个独立的系数 ...

大哥你的意思是不是说，我最后得到的值不是方差的值，而是方差*（n-1）的值？

xianshasaman

对，varOut里面就是  方差*（n-1）