LFOSSA 技能焕新季即将启动： 7 折福利超前锁定，仅限前 30 名！

描述统计

描述统计是研究

1. 如何取得反映客观现象的数据(数据的收集)    
2. 通过图表形式对数据进行加工处理和可视化。
3. 通过概括与分析得出反应客观现象的规律性数量特征。

数据的可靠性和有效性

数据是否可靠(reliable)和有效(valid)?

1. 可靠性(reliability):多次测量得到的数据是否一致。
2. 有效性(validity):实际测量的对象=认为/希望测量的对象

分类变量特征和可视化

• 无序分类变量

性别(名目；=,≠):男，女

观测12个新生儿的性别(n=12),结果为：女，男，女，女，男，男，男，男，女，男，男，女

频率表(frequency table)

条形图(bar plot)

如果我们将频数换成频率，那么左边的刻度值就是百分比，女性的上限频率值就是41.7%，男性的上限频率值就是58.3%。

集中趋势(central tendency):一组观测值向其中心集中的倾向和程度，对于无序分类来说只有一种

众数(mode):一组观测值中出现次数最多的数，对于上面的男女分类来说，由于男性的数量更多，故它的众数就是男。众数主要体现在大家关注的大多数，比如考试完对选择题的答案，如果你选择的跟大多数人一致，那么你可能就会认为你选对了。但一组观测值中可能存在多个众数，也可能不存在众数，比如颜色：赤1，橙1，黄1，绿1，青1，蓝1，紫1，此时不存在众数；赤2，橙6，黄1，绿10，青3，蓝10，紫4，此时存在多个众数。

• 有序分类变量

有序分类变量可以进行大小的比较，但是不能进行加减法。

教育程度(次序；=,≠，>，<)

小学(1)，初中(2)，高中(3)，本科(4)，研究生(5)

观测19个人的教育程度(n=19)，结果为：3，3，4，1，5，4，2，1，5，4，4，4，5，3，2，1，4，5，5

频率表

条形图

集中趋势：众数  ->  本科

集中趋势：中位数(median)：对于有限的数集，把所有的观测值按大小排序后，位于正中间的观测值即为中位数/中值。

我们将上面教育程度的数值进行排序：1，1，1，2，2，3，3，3，4，4，4，4，4，4，5，5，5，5，5

所以这一组数的中位数就是位于第10位的4。以上是19个数，它是一个奇数，所以只有一个中位数，如果为偶数，如：1，1，1，2，2，3，3，3，4，4，4，4，4，4，5，5，5，5，5，5

那么我们会将第10位的4和第11位的4加和除以2，即\({4+4\over 2}=4\)

数值变量特征和可视化

• 等距变量

温度(=,≠,>,<,+,-)

5月份前两周的温度：19，22，21，17，13，19，18，17，17，21，21，21，19，20   (n=14)

频率表

数值变量与分类变量频率表的区别，数值变量等距，可以分割小区间，令Δ=1

频率直方图(histogram)

频率直方图的横坐标是温度，纵坐标为频率/Δ，注意这里不是频率，这么做的目的是使蓝色覆盖的区域面积为1，方便进行归一化。

若Δ=5，有

频率直方图

集中趋势：众数、中位数，我们先将这些温度进行排序

13，17，17，17，18，19，19，19，20，21，21，21，21，22

众数=21

中位数=\({19+19\over 2}=19\)

对于数值变量，由于可以做加减法，于是有了第三个集中趋势指标——均值(mean)：在一组数据中，所有数据之和再除以这组数据的个数，所得即为这组数据的均值。

=\({\sum_{i=1}^nX_i\over n}\)

上面的这组温度值的均值就为

\({19+22+21+17+13+19+18+17+17+21+21+21+19+20\over 14}={265\over 14}=18.93\)

• 离散趋势(tendency of dispersion)：观测值偏离其中心的趋势。

1、级差/全距(Range)：最大值减去最小值，用于简单描述数据的范围大小。

上面的这组温度值的级差就为

级差=22-13=9

级差的缺点是非常容易受到极端值的影响，比如在这组温度值中的13在某种程度上就是一个极端值。

2、分位数/分位点(quantile)：把数据n等分的分割点

之前我们知道这组温度的中位数 =\({19+19\over 2}=19\)，中位数把数据分成了数目相等的两部分(排序后的)，是二分位数/点。

分位数常见的有四分位数(quartile)

我们将中位数的左右两边继续两等分，就有了17和21两个四等分点，这三个四等分点分别称为\(Q_1\)(25%分位点)、\(Q_2\)(50%分位点)、\(Q_3\)(75%分位点)。\(Q_1\)和\(Q_3\)的距离称为四分位距(interquartile range,IQR)，\(IQR=Q_3-Q_1\)，因此在这组温度值中IQR=21-17=4。

如果中位数的左右两边的数据个数为偶数个，得到\(Q_1\),\(Q_3\)的方式与之前得到中位数的方式相同。

3、箱线图(box plot)

箱线图不是完全对称的，在上图中，最小值距离中位数的距离比最大值距离中位数的距离要长，说明在中位数之下的这些数分布的范围更广，中位数之上的数分布的范围更加集中。箱线图提供了更多在离散趋势中的度量，一旦有了箱线图就可以对数据的分布有了一定的了解。箱线图也可以帮助我们发现极端值。

• 等比变量

鸢尾花花瓣的长度(=,≠,>,<,+,-,*,/)

鸢尾花花瓣能够进行乘除运算的原因是鸢尾花花瓣具有真0值，长度0代表鸢尾花花瓣不存在。虽然之前的等距变量的温度也有0度，但是不代表这个温度不存在，故等距变量没有乘除法。

假设我们观测了50朵鸢尾花的花瓣，直接给出频数表(n=50)

频率直方图(Δ=0.1)

众数=5.0和5.1

级差=5.8-4.3=1.5

均值=\({4.3*1+4.4*3+4.5*1+4.6*4+4.7*2+4.8*5+4.9*4+5.0*8+5.1*8+5.2*3+5.3*1+5.4*5+5.5*2+5.7*2+5.8*1\over 50}=5.0\)

中位数=5.0

箱线图

IQR=5.2-4.8=0.4

离散趋势：方差(variance)和标准差(standard deviation)

方差：每一个观测值与均值之间的差异的平方和的平均数

\(σ^2={\sum_{i=1}^n(X_i-mean(X))^2\over n}\)

这组鸢尾花的方差=\({(4.3-5.0)^2+3*(4.4-5.0)^2+(4.5-5.0)^2+4*(4.6-5.0)^2+2*(4.7-5.0)^2+5*(4.8-5.0)^2+4*(4.9-5.0)^2+8*(5.0-5.0)^2+8*(5.1-5.0)^2+3*(5.2-5.0)^2+(5.3-5.0)^2+5*(5.4-5.0)^2+2*(5.5-5.0)^2+2*(5.7-5.0)^2+(5.8-5.0)^2\over 50}\)

=0.124

标准差：方差直接开方

\(σ=\sqrt{σ^2}=\sqrt{{\sum_{i=1}^n(X_i-mean(X))^2\over n}}\)

这组鸢尾花的标准差=\(\sqrt{0.124}=0.352\)

标准差与原观测值具有相同的单位

由此，表现集中趋势的主要有：众数、中位数、均值；表现离散趋势的主要有：级差、分位数、方差、标准差。

分布的形状

在鸢尾花的频率直方图

中为我们提供了很多的信息。通过该图，我们知道了出现次数最多的鸢尾花花瓣的长度为5厘米，它的变化范围为4.3厘米到5.8厘米之间。

箱线图也为我们提供了有用的信息

紫色箱子上方的虚线存在着25%的数据，箱子本身存在着50%的数据，箱子下方的虚线存在着25%的数据。

当我们看到某一个数据的分布的时候，才能够判断它是不是跟我们预期中的一样，是否符合某个已知的概率分布，从而决定我们使用哪种统计方法对其进行分析。

• 偏度(skewness)

在上图的第一个图中，我们可以看到在它的左边有一个长尾，如果我们把该图看成是人类的预期寿命的话，那么我们知道预期寿命短的人是很少的，大部分人的寿命都集中在某一个峰值。通过该图，我们可以得到众数，绿线为中位数，红线为均值。在左侧的长尾中，均值会被往左拉，均值受到了特别小的数值影响会小于中位数，中位数是排序后出现在中间位置的数，它不受左边长尾小数值的影响。对于这种分布，我们称之为左偏(均值<中位数≤众数)。

在上图的第二个图中，它是比较对称的，现实生活中，有很多观测数据是比较符合这种分布的，比如人类的身高和体重。在该图中，红线和绿线是重叠的，对于这种分布，我们称之为对称(均值=中位数=众数)。

在上图的第三个图中，它与第一个图的分布是相反的，它的右侧有一个长尾，如果我们把该图看成是某公司员工薪水的话，那么我们知道大部分人的收入较低，高薪的人是较少的。在右侧的长尾中，均值会被往右拉，受到了高薪数值的影响会大于中位数，中位数也是薪资排序后出现在中间位置的数，不受高薪数值的影响。对于这种分布，我们称之为右偏(均值>中位数≥众数)。

• 形态(modality)

形态是用来描述分布有多少个峰。上图中的第一个图只有一个峰，这种分布我们称之为单峰(unimodal)，我们的身高和体重就符合这种分布；上图中的第二个图有两个峰，这种分布我们称之为双峰(bimodal)，假设我们同时收集儿童和成年人的身高，那么儿童的身高会集中在某个比较低的数值，成年人的身高会集中在某个比较高的数值，如果我们的数据出现了双峰，则我们需要分析出现的原因。上图中的第三个图中有三个峰，这种有三个峰或者多于三个峰的分布，我们称之为多峰(multimodal)，多峰分布其实并不常见，当我们看到多峰的时候需要分析这是反应真实客观的分布情况还是我们数据采集出了问题。

• 峰度(kurtosis)

该图是在频率直方图的基础上，将数值进行了连线，形成的一个数值走向曲线。通过该图，我们可以看到这组数据峰尖、尾平，数据向中心聚拢的程度是非常高的。

与中心聚拢程度高相反的是这种扁平的，数据向中心聚拢程度低的，数据的分布是相对分散的。

在该图中，我们可以看到它的峰没有那么尖锐，也没有那么扁平，大部分数据处于峰的附近，说明数据向中心聚拢的程度是比较高的。

变量间的关系

• 两个分类变量的关系

泰坦尼克号的325名乘客：幸存(是/否)；年龄(儿童/成人)

年龄与幸存是否有关系？

1、关联表(contingency table)

相对频率表

从以上的数据显示它们似乎有关系，又似乎没有关系，我们需要经过一个统计推断的过程才能得出正确的结论。这里需要注意的是我们说两个变量有依赖关系不等于说它们有因果关系。

2、分段条形图

上图中，横坐标是年龄这个变量，颜色是幸存与否这个变量，绿色表示幸存，紫色表示未幸存。纵坐标表示人数。通过这个图形，我们可以看到，儿童的总数与成年人的总数是有巨大差异的。通过这个图形，我们也可以初步判断年龄与幸存是有关系的，但并不能作为结论。

相对频率分段条形图

这里的纵坐标不再是人数而是百分比，这样对于儿童和成人，它们的群体都是100%，所以是一样高的。同样，绿色代表幸存，紫色代表未幸存。

分段条形图和相对频率分段条形图有略微的不同，从分段条形图中，我们可以看出儿童和成人在总数上的差异；而相对频率分段条形图更直观呈现了幸存与否在儿童和成人之间的占比。

• 两个数值变量的关系

工资与入职时间是否有关？什么关系？关系强弱？房价与到学校的距离是否有关？什么关系？关系强弱？

散点图(scatter plot)：方向、形状、强度、极端值

在上图中，横坐标代表入职时间，纵坐标代表工资，我们可以看到有些点的入职时间很短，但是工资却很高，而有些点的入职时间很长，但是工资却不是很高，这意味着两个变量有关系的可能性很小。

在上图中，我们可以看到所有的点都在一条直线上。这意味着两个变量是完美的线性相关，随着一个变量的增加，另一个变量一定会增加。这个增加代表的是关系的方向；由于是一条线，是线性关系，代表的是两个数值变量关系的形状；由于排列在一条直线上，说明两个变量的关系强度是非常强的。

在这张图中，虽然它的分布也比较分散，但是相对于第一张图，它还是有着相对明显的趋势的，整体的趋势是随着第一个变量的增加，第二个变量也增加，并且两者是线性关系，但是它的强度显然没有第二张图中两者的关系强。其实这种情况才往往是我们在现实中看到的情况。因为我们采集到的数据往往是有误差和噪音的，几乎不可能得到一条类似于第二张图的情况。

这张图与第三张图是类似的，只不过两个变量之间是负相关。

在第二、第三、第四张图中，两个变量之间都是线性相关的，但在这张图中，这些点形成了一个曲线的分布，这意味着两个数值变量之间的关系不是简单的单调递增或者单调递减，而是有一个拐点，在某个点之前是递减，在某个点之后是递增，并且不是线性的递增，而是非线性的递增。

在这张图中，我们可以看到有两个点非常远离了左下角的散点，则我们可以称这两个点为极端值。如果在我们的图形中出现了这样的点，一定要引起我们的注意，极端值的存在会影响结论的真实性，比如说本来左下角这一组散点就像第一张图一样，两个变量是无关的，但是加上了右上角这个点之后，我们就会得到一个假的关系，随着一个变量的增大，另一个变量也会增大，这个虚假的关系就是由这一个极端值导致的。

• 一个数值变量和一个分类变量的关系

工资与性别的关系，房价与学区房的关系

并排箱图(side-by-side box plot)

上图中横坐标就是性别的分类变量，纵坐标就是工资的数值变量，对于这张图来看，男性的工资的中位数是远远高于女性工资的中位数的，由该图我们可以初步得出工资水平是和性别有关系的。

极端值与缺失值

• 极端值/异常值(outliers)

1. 在一组数据中，小于\(Q_1-1.5IQR\)或者大于\(Q_3+1.5IQR\)的数据是疑似极端值。
2. 在一组数据中，小于\(Q_1-3IQR\)或者大于\(Q_3+3IQR\)的数据是极端值。

修正箱线图

极端值的产生可能源于

1. 数据的测量、记录或输入时的错误；
2. 数据来自不同的总体(例如：病人 vs 健康人)，这里就好像之前的双峰时候的儿童与成年人的身高数据；
3. 数据是正确的，但它只体现小概率事件；

极端值可能产生的影响

某公司员工的收入水平，1.2,1.3,1.4,1.5,1.6,1.6,1.8,2.0,2.2,15

当我们去计算均值的时候，如果包含了最后的15w，那么此时均值为2.96w；如果剔除这个15w，那么此时的均值为1.62w。此时我们可以看到这个15w远远拉高了均值，会给人一种错误。我们看看中位数，无论是否有这个15w，中位数都是1.6w，它非常接近不带极端值的均值。我们看看众数为1.6w，那么我们就知道，均值受异常值的影响非常大，而中位数和众数几乎不受异常值的影响。当我们在计算均值之前一定要处理极端值的情况。这些都是集中趋势的指标。

现在我们来看看离散趋势的指标，级差：如果带15w这个数据，那么级差为13.8w，如果去掉15w，那么级差为1w。由此可见级差也是对极端值非常的敏感，在计算级差值之前也要处理掉极端值。再来看标准差，如果带15w的话，标准差为4.24，不带15w的话，标准差为0.33，由此可见标准差也是受极端值影响非常严重的。受极端值影响小的离散趋势的指标为IQR，这里的IQR为，带15w的IQR=2.0-1.4=0.6，不带15w的IQR=1.9-1.35=0.55。

如何处理极端值

1. 如果是由于测量或记录的错误，或其他明显的原因造成的，直接丢弃即可。
2. 如果极端值出现的原因无法解释，那么，丢弃或保留极端值则需要具体问题具体分析；尽量选用受极端值影响小的指标，如中位数、众数、IQR。
3. 可以通过对比保留极端值和丢弃极端值对结果的影响，来判断结果是否受到极端值的影响。

• 缺失值

如何处理缺失值

1. 如果含有缺失值的观测记录很少，而数据量很大，可以把含有缺失值的观测记录丢弃。如果当我们要观测1000个人的身高，但其中有10个人的记录丢失，由于从1000个人和从990个人的身上得到的结论区别不大，此时我们可以丢弃这10个人的缺失值。
2. 如果含有缺失值的观测记录很多，需要分析原因，看是否能够把缺失的记录补全。
3. 如果含有缺失值的观测记录较少，可以使用均值/中位数/众数/最大值来进行替代。

描述统计编程实现

频数

from collections import Counter

if __name__ == '__main__':

    data = [2, 2, 2, 2, 1, 1, 1, 3, 3]  # 样本
    counter = Counter(data)
    print(counter.most_common())  # 频数

运行结果

[(2, 4), (1, 3), (3, 2)]

从结果我们可以看出，在样本中，2有4个，1有3个，3有2个。

numpy实现

import numpy as np

if __name__ == '__main__':

    data = np.array([2, 2, 2, 2, 1, 1, 1, 3, 3])
    unique_elements, counts = np.unique(data, return_counts=True)
    print(dict(zip(unique_elements, counts)))

运行结果

{1: 3, 2: 4, 3: 2}

pytorch实现

import torch

if __name__ == '__main__':

    data = torch.tensor([2, 2, 2, 2, 1, 1, 1, 3, 3])
    unique_values, counts = torch.unique(data, return_counts=True)
    print(dict(zip(unique_values, counts)))

运行结果

{tensor(1): tensor(3), tensor(2): tensor(4), tensor(3): tensor(2)}

频率

from collections import Counter

if __name__ == '__main__':

    data = [2, 2, 2, 2, 1, 1, 1, 3, 3]  # 样本
    counter = Counter(data)
    ret = []
    for point in counter.most_common():
        ret.append((point[0], point[1] / len(data)))  # 频率
    print(ret)

运行结果

[(2, 0.4444444444444444), (1, 0.3333333333333333), (3, 0.2222222222222222)]

从结果我们可以看出，在样本中，2的占比为44.4%，1的占比为33.3%，3的占比为22.2%。

numpy实现

import numpy as np

if __name__ == '__main__':

    data = np.array([2, 2, 2, 2, 1, 1, 1, 3, 3])
    unique_elements, counts = np.unique(data, return_counts=True)
    print(dict(zip(unique_elements, counts / len(data))))

运行结果

{1: 0.3333333333333333, 2: 0.4444444444444444, 3: 0.2222222222222222}

pytorch实现

import torch

if __name__ == '__main__':

    data = torch.tensor([2, 2, 2, 2, 1, 1, 1, 3, 3])
    unique_values, counts = torch.unique(data, return_counts=True)
    print(dict(zip(unique_values, counts / len(data))))

运行结果

{tensor(1): tensor(0.3333), tensor(2): tensor(0.4444), tensor(3): tensor(0.2222)}

集中趋势

• 众数

from collections import Counter

if __name__ == '__main__':

    data = [2, 2, 2, 2, 1, 1, 1, 1, 3, 3]  # 样本
    counter = Counter(data)
    ret = []
    if counter.most_common()[0][1] == 1:
        print('众数不存在')
    else:
        count = counter.most_common()[0][1]  # 获取样本中最多的元素的个数
        for point in counter.most_common():
            if point[1] == count:  # 查看有多少元素拥有相同最多的个数
                ret.append(point[0])
            else:
                break
        print(ret)  # 打印众数

运行结果

[2, 1]

在这里，我在样本中比之前多加了一个1，从结果我们可以看出，2和1的个数最多为4个，所以2和1都为众数。

numpy实现

import numpy as np

if __name__ == '__main__':

    data = np.array([2, 2, 2, 2, 1, 1, 1, 1, 3, 3])
    ret = []
    unique_elements, indices, counts = np.unique(data, return_index=True, return_counts=True)
    max_count = np.max(counts)
    mode = indices[counts == max_count]
    for i in mode:
        ret.append(data[i])
    print(ret)

运行结果

[1, 2]

pytorch实现

在Pytorch中，如果样本只有一个众数，可以使用torch.mode()来获取

import torch

if __name__ == '__main__':

    data = torch.tensor([2, 2, 2, 2, 1, 1, 1, 3, 3])
    print(torch.mode(data))

运行结果

torch.return_types.mode(
values=tensor(2),
indices=tensor(1))

但是如果有多个众数的话，得自己实现

import torch

if __name__ == '__main__':

    data = torch.tensor([2, 2, 2, 2, 1, 1, 1, 1, 3, 3])
    ret = []
    unique_values, counts = torch.unique(data, return_counts=True)
    max_count = torch.max(counts)
    for i, count in enumerate(counts):
        if count == max_count:
            ret.append(unique_values[i])
    print(ret)

运行结果

[tensor(1), tensor(2)]

•  中位数

from collections import Counter

if __name__ == '__main__':

    data = [2, 2, 2, 2, 1, 1, 1, 1, 3, 3]  # 样本
    sorted_data = sorted(data)  # 对样本进行排序
    n = len(data)
    if n % 2 == 1:  # 如果样本有奇数个
        print(sorted_data[n // 2])  # 中位数
    else:  # 样本有偶数个的中位数
        print((sorted_data[n // 2 - 1] + sorted_data[n // 2]) / 2)

运行结果

2.0

numpy实现

import numpy as np

if __name__ == '__main__':

    data = np.array([2, 2, 2, 2, 1, 1, 1, 1, 3, 3])
    print(np.median(data))

运行结果

2.0

pytorch实现

import torch

if __name__ == '__main__':

    data = torch.tensor([2, 2, 2, 2, 1, 1, 1, 1, 3, 3])
    print(torch.median(data))

运行结果

tensor(2)

• 均值

from collections import Counter

if __name__ == '__main__':

    data = [2, 2, 2, 2, 1, 1, 1, 1, 3, 3]  # 样本
    print(sum(data) / len(data))  # 均值

运行结果

1.8

numpy实现

import numpy as np

if __name__ == '__main__':

    data = np.array([2, 2, 2, 2, 1, 1, 1, 1, 3, 3])
    print(np.mean(data))

运行结果

1.8

pytorch实现

import torch

if __name__ == '__main__':

    data = torch.tensor([2, 2, 2, 2, 1, 1, 1, 1, 3, 3], dtype=torch.float32)
    print(torch.mean(data))

运行结果

tensor(1.8000)

离散趋势

• 级差

from collections import Counter

if __name__ == '__main__':

    data = [2, 2, 2, 2, 1, 1, 1, 1, 3, 3]  # 样本
    print(max(data) - min(data)) # 级差

运行结果

2

numpy实现

import numpy as np

if __name__ == '__main__':

    data = np.array([2, 2, 2, 2, 1, 1, 1, 1, 3, 3])
    print(np.ptp(data))

运行结果

2

pytorch实现

import torch

if __name__ == '__main__':

    data = torch.tensor([2, 2, 2, 2, 1, 1, 1, 1, 3, 3])
    print(torch.max(data) - torch.min(data))

运行结果

tensor(2)

• 四分位数

from collections import Counter

def median(data):
    n = len(data)
    if n % 2 == 1:
        return data[n // 2]
    else:
        return (data[n // 2 - 1] + data[n // 2]) / 2

if __name__ == '__main__':

    data = [2, 2, 2, 2, 1, 1, 1, 1, 3, 3]  # 样本
    n = len(data)
    if n >= 4:
        sorted_data = sorted(data)  # 对样本进行排序
        q2 = median(sorted_data)  # 获取中位数，即50%分位点
        if n % 2 == 1:   # 如果样本有奇数个
            q1 = median(sorted_data[:n // 2])  # 获取25%分位点
            q3 = median(sorted_data[n // 2 + 1:])  # 获取75%分位点
        else:  # 如果样本有偶数个
            q1 = median(sorted_data[:n // 2])  # 获取25%分位点
            q3 = median(sorted_data[n // 2:])  # 获取75%分位点
        print(q1, q2, q3)

运行结果

1 2.0 2

numpy实现

import numpy as np

if __name__ == '__main__':

    data = np.array([2, 2, 2, 2, 1, 1, 1, 1, 3, 3])
    q1 = np.percentile(data, 25)
    q2 = np.percentile(data, 50)
    q3 = np.percentile(data, 75)
    print(q1, q2, q3)

运行结果

1.0 2.0 2.0

pytorch实现

import torch

if __name__ == '__main__':

    data = torch.tensor([2, 2, 2, 2, 1, 1, 1, 1, 3, 3], dtype=torch.float32)
    q1 = torch.quantile(data, 0.25)
    q2 = torch.quantile(data, 0.5)
    q3 = torch.quantile(data, 0.75)
    print(q1, q2, q3)

运行结果

tensor(1.) tensor(2.) tensor(2.)

• 方差

from collections import Counter

if __name__ == '__main__':

    data = [2, 2, 2, 2, 1, 1, 1, 1, 3, 3]  # 样本
    n = len(data)
    if n > 1:
        mean_value = sum(data) / n  # 均值
        var = sum((e - mean_value)**2 for e in data) / n  # 方差
        print(var)

运行结果

0.5599999999999999

numpy实现

import numpy as np

if __name__ == '__main__':

    data = np.array([2, 2, 2, 2, 1, 1, 1, 1, 3, 3])
    print(np.var(data))

运行结果

0.5599999999999999

pytorch实现

import torch

if __name__ == '__main__':

    data = torch.tensor([2, 2, 2, 2, 1, 1, 1, 1, 3, 3], dtype=torch.float32)
    print(torch.var(data, unbiased=False))

运行结果

tensor(0.5600)

• 标准差

from collections import Counter
import math

if __name__ == '__main__':

    data = [2, 2, 2, 2, 1, 1, 1, 1, 3, 3]  # 样本
    n = len(data)
    if n > 1:
        mean_value = sum(data) / n  # 均值
        var = sum((e - mean_value)**2 for e in data) / n  # 方差
        std = math.sqrt(var)  # 标准差
        print(std)

运行结果

0.7483314773547882

numpy实现

import numpy as np

if __name__ == '__main__':

    data = np.array([2, 2, 2, 2, 1, 1, 1, 1, 3, 3])
    print(np.std(data))

运行结果

0.7483314773547882

pytorch实现

import torch

if __name__ == '__main__':

    data = torch.tensor([2, 2, 2, 2, 1, 1, 1, 1, 3, 3], dtype=torch.float32)
    print(torch.std(data, unbiased=False))

运行结果

tensor(0.7483)

作图

• 散点图

import matplotlib.pyplot as plt
import random

if __name__ == '__main__':

    random.seed(666)
    x = [random.randint(0, 100) for _ in range(100)]  # 获取100个0～100的随机数
    y = [random.randint(0, 100) for _ in range(100)]
    plt.scatter(x, y)  # 散点图
    plt.show()

运行结果

• 折线图

import matplotlib.pyplot as plt
import random

if __name__ == '__main__':

    random.seed(666)
    x = [random.randint(0, 100) for _ in range(100)]  # 获取100个0～100的随机数
    plt.plot([i for i in range(100)], x)  # 折线图
    plt.show()

运行结果

• 条形图

import matplotlib.pyplot as plt
from collections import Counter

if __name__ == '__main__':

    data = [3, 3, 4, 1, 5, 4, 2, 1, 5, 4, 4, 4, 5, 3, 2, 1, 4, 5, 5]
    counter = Counter(data)
    x = [point[0] for point in counter.most_common()]  # 获取元素
    y = [point[1] for point in counter.most_common()]  # 获取频数
    plt.bar(x, y)  # 条形图
    plt.show()

运行结果

这个是代表频数的条形图，如果要绘制频率的条形图，只需要修改程序为提取频率的数据即可

import matplotlib.pyplot as plt
from collections import Counter

if __name__ == '__main__':

    data = [3, 3, 4, 1, 5, 4, 2, 1, 5, 4, 4, 4, 5, 3, 2, 1, 4, 5, 5]
    counter = Counter(data)
    x = [point[0] for point in counter.most_common()]  # 获取元素
    y = [point[1] / len(data) for point in counter.most_common()]  # 获取频率
    plt.bar(x, y)  # 条形图
    plt.show()

运行结果

• 直方图

import matplotlib.pyplot as plt
import random

if __name__ == '__main__':

    random.seed(666)
    data = [random.randint(1, 100) for _ in range(1000)]  # 获取1000个1～100的随机数
    plt.hist(data, rwidth=0.8, bins=10, density=True)  # 频率直方图
    plt.show()

运行结果

• 箱线图

import matplotlib.pyplot as plt
import random

if __name__ == '__main__':

    random.seed(666)
    data = [random.randint(1, 100) for _ in range(1000)]  # 获取1000个1～100的随机数
    plt.boxplot(data)  # 箱线图
    plt.show()

运行结果

这是一个没有极端值的数据，现在我们来添加极端值

import matplotlib.pyplot as plt
import random

if __name__ == '__main__':

    random.seed(666)
    data = [random.randint(1, 100) for _ in range(1000)]  # 获取1000个1～100的随机数
    data.append(200)
    data.append(-200)
    plt.boxplot(data)  # 箱线图
    plt.show()

运行结果

并排箱图

import matplotlib.pyplot as plt
import random

if __name__ == '__main__':

    random.seed(666)
    data1 = [random.randint(66, 166) for _ in range(200)]  # 获取200个66～166的随机数
    data2 = [random.randint(60, 120) for _ in range(200)]  # 获取200个60～120的随机数
    plt.boxplot([data1, data2])  # 并排箱图
    plt.show()

运行结果

样本和抽样分布

之前都属于描述统计的范畴，这里我们开始进入统计推断的阶段，统计推断又包括——参数估计和假设检验。

总体(全样本)与样本(抽样样本)

总体又叫做全样本，如果这个对象是人的话，那么它包含了这个世界上所有的人，但往往我们作统计分析的时候不可能拿到全世界所有人的数据，故我们会在某些地域采集部分的人的数据，称为样本，我们可能是分批采集的不同地域的样本，如上图右边所示。这里的每一份样本都可以向我们提供一些总体的信息。参数估计就是估计的一些总体的参数，不同的样本提供的信息是有一定差别的。我们如何使用样本的信息去得到总体的信息，解决这个问题的过程就叫做参数估计。我们得到的估计是否正确，就是假设检验要解决的问题。

1. 总体：试验的全部可能的观察值。
2. 个体：每一个观察值。
3. 总体的容量：总体中所包含的个体的个数。
4. 一个总体对应于一个随机变量X。
5. 对总体的研究就是对一个随机变量X的研究，X的分布函数和数字特征就称为总体的分布函数和数字特征。
6. 在实际中，总体的分布一般是未知的，或只知道它具有某种形式而其中包含着未知参数。
7. 从总体中抽取一部分个体，根据获得的数据来对总体分布做出推断。被抽出的部分个体叫做总体的一个样本。
8. 样本是进行统计推断的依据。
9. 从总体抽取一个个体，就是对总体X进行一次观察并记录其结果。
10. 在相同的条件下对总体X进行n次重复的、独立的观察，将n次观察结果按试验的次序记为\(X_1,X_2,...,X_n\)，那么这组数据就称为来自总体X的一个简单随机样本，n称为样本容量。\(X_1,X_2,...,X_n\)是独立且与X同分布的，这意味着这里的每一个\(X_i\)(i=1,...,n)就是一个随机变量，它可能的取值和取值出现的概率都跟总体X是一致的。其原因为\(X_i\)代表的是X中的某一个个体，在我们真正从这个个体中得到观察值之前，它可能的取值跟总体可能的取值范围是一致的，并且它所对应的观察值可能出现的概率也跟总体的情况是一致的，因此我们认为\(X_1,X_2,...,X_n\)它们每一个都是随机变量，并且与总体具有相同的分布。
11. 当n次观察一经完成，我们就得到一组实数\(x_1.x_2,...,x_n\)，它们是\(X_1,X_2,...,X_n\)的观察值，称为样本值。

抽样分布

我们往往不是直接使用样本本身，而是针对不同的问题构造样本的适当函数，利用这些样本的函数进行统计推断。这个样本的函数称为统计量。

• 统计量

\(X_1,X_2,...,X_n\)是来自总体X的一个样本，\(g(X_1,X_2,...,X_n)\)是\(X_1,X_2,...,X_n\)的函数，若g中不含未知参数，则称\(g(X_1,X_2,...,X_n)\)是一个统计量。

常用统计量

\(X_1,X_2,...,X_n\)是来自总体X的一个样本，\(x_1.x_2,...,x_n\)是这一样本的样本值。

1. 样本均值：\(g(X_1,X_2,...,X_n)\)= =\({\sum_{i=1}^nX_i\over n}\)，其观察值为\(g(x_1,x_2,...,x_n)\)==\({\sum_{i=1}^nx_i\over n}\)
2. 样本方差：\(S^2={\sum_{i=1}^n(X_i-mean(X))^2\over n-1}\)，其观察值为\(S^2={\sum_{i=1}^n(x_i-mean(x))^2\over n-1}\)，这里需要说明的是由于我们使用的是样本而并不是总体来求方差，除以n-1是为了得到无偏的估计，而除以n的话得到的就是有偏的估计。
3. 样本标准差：\(S=\sqrt{S^2}\)

随着样本统计量的观察值被获取，它们也会形成一个分布，这个分布我们就叫做统计量的分布，也就是我们所说的抽样分布。

常用统计量的分布

• 来自正态分布总体的几个常用统计量分布

1、卡方分布(chi-square distribution)

总体服从标准正态分布N(0,1)，\(X_1,X_2,...,X_n\)是来自总体N(0,1)的一个样本，我们需要的统计量卡方分布为

\(χ^2=X_1^2+X_2^2+...+X_n^2\)

服从自由度为n的卡方分布记为\(χ^2\)~\(χ^2(n)\)，自由度(df:degree of freedom)指的是独立变量的个数。

卡方分布的概率密度函数

\(f_n(x)={x^{{n\over 2}-1}e^{-{x\over 2}}\over 2^{n\over 2}Gamma({n\over 2})\sqrt{π}}\)

其中\(Gamma({n\over 2})\)是伽马函数。

卡方分布的概率密度函数图

当自由度为1的时候，也就是一次只取出一个个体，我们可以看到当取值比较小接近0时的概率是最大的。

当自由度为3的时候，概率密度最大的值比n=1的时候要向右偏移了一些。

当自由度为6的时候，卡方分布图更加的胖了一些，并且概率密度最大的值更加向右偏移，到达了5附近。

当自由度为12的时候，卡方分布图更加的胖了一些，并且概率密度最大的值更加向右偏移，到达了10附近。

当自由度为24的时候，在我们的可见范围内，概率密度最大的值到达了20。

由此我们可以知道，同样是基于样本构建出来的卡方这个统计量，它的分布的样子取决于样本容量即自由度。

卡方分布的性质

1. 可加性：\(χ_1^2\)~\(χ^2(n_1)\)，\(χ_2^2\)~\(χ^2(n_2)\)，并且\(χ_1^2\),\(χ_2^2\)相互独立，则有\(χ_1^2+χ_2^2\)~\(χ^2(n_1+n_2)\)
2. 期望(均值)和方差：\(E(χ^2)=n\)，\(Var(χ^2)=2n\)

用Python绘画卡方分布概率密度函数图

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma


# 定义卡方分布的概率密度函数
def chi_squared_pdf(x, n):
    return (x**(n / 2 - 1) * np.exp(-x / 2)) / (gamma(n / 2) * 2**(n / 2) * np.sqrt(np.pi))


if __name__ == '__main__':

    # 设置要绘制的卡方分布的自由度n
    n = 5

    # 创建1000个从0.01到30的等差数列
    x = np.linspace(0.01, 30, 1000)

    # 计算卡方分布的概率密度函数
    pdf_values = chi_squared_pdf(x, n)

    # 绘制卡方分布概率密度函数曲线图
    plt.plot(x, pdf_values, label='n = {}'.format(n))

    # 设置图表标题和轴标签
    plt.title('Chi-Squared Probability Density Function for n = {}'.format(n))
    plt.xlabel('x')
    plt.ylabel('Probability Density')

    # 显示图例
    plt.legend()

    # 显示图表
    plt.show()

运行结果

2、t分布(t disttribution)

总体服从标准正态分布N(0,1)。即样本统计量X~N(0,1)，样本统计量Y服从自由度为n的卡方分布，即Y~\(χ^2(n)\)，X、Y相互独立，我们需要的统计量t分布为

\(t={X\over \sqrt{{Y\over n}}}\)

服从自由度为n的t分布记为t~t(n)，t分布又称为Student分布。

这里需要注意的是，这里的X、Y并不是指的随机变量，而是样本统计量，X是由\(X_1,X_2,...,X_n\)得到的均值，Y是由\(X_1^2+X_2^2+...+X_n^2\)得到的方差。

t分布的概率密度函数

\(f_t(x)={Gamma({n+1\over 2})\over \sqrt{nπ}Gamma({n\over 2})}(1+{x^2\over n})^{-{n+1\over 2}}\)

t分布的概率密度函数图

当自由度为1的时候，上图中的黑线是标准正态分布图，红线就是自由度为1的t分布的图，我们可以看到t分布比较像标准正态分布，但是它的尾部(-4到-2，2到4的部分)比标准正态分布要宽，它的尖部比标准正态分布要窄。

当自由度为3的时候，上图中的黑线依然是标准正态分布图，绿线是自由度为3的t分布的图，我们可以看到t分布的图更加接近于标准正态分布，与n=1的t分布相比，绿线与黑线更加贴近。

当自由度为10的时候，蓝线所代表的t分布的图与黑线的标准正态分布图已经极其的接近了，只在0附近、-4到-2、2到4留了一点点缝隙，但这个缝隙已经远远小于n=1和n=3时到缝隙了。

当自由度为20的时候，黄线所代表的t分布的图已经几乎和黑线的标准正态分布图重合了。

当自由度为45的时候，粉线所代表的t分布的图几乎已经和黑线的标准正态分布图完全重合了。

用Python绘画t分布概率密度函数图

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import t

if __name__ == '__main__':

    # 定义t分布的参数
    n = 10  # 自由度
    loc = 0  # 位置参数
    scale = 1  # 比例参数

    # 创建从200个从-5到5的等差数列
    x = np.linspace(-5, 5, 200)

    # 计算t分布的概率密度值
    t_pdf = t.pdf(x, n, loc, scale)

    # 绘图
    plt.plot(x, t_pdf, 'b-', label='t-distrib, $n={}$, loc={}, scale={}'.format(n, loc, scale))

    # 设置轴标签和标题
    plt.title('t-distribution PDF')
    plt.xlabel('X')
    plt.ylabel('Probability Density')

    # 显示图例
    plt.legend()

    # 显示图形
    plt.show()

运行结果

3、F分布(F distribution)

总体服从标准正态分布 N (0,1)，即样本统计量 X~N (0,1)，样本统计量U服从自由度为\(n_1\)的卡方分布，即U~\(χ^2(n_1)\)，样本统计量V服从自由度为\(n_2\)的卡方分布，即V~\(χ^2(n_2)\)，U、V相互独立，我们需要的统计量F分布为

\(F={{U\over n_1}\over {V\over n_2}}\)

服从自由度为\((n_1,n_2)\)的F分布，记为F~\(F(n_1,n_2)\)

F分布的概率密度函数

F分布的概率密度函数图

上图是固定了\(n_2\)为40，变换\(n_1\)的F分布的概率密度函数图。我们可以看到在\(n_2\)固定的情况下，随着\(n_1\)的减小，F分布的图形变的越来越矮，越来越胖，并且其概率密度最大值所对应的F值越来越向左偏移。

上图是固定了\(n_1\)为20，变换\(n_2\)的F分布的概率密度函数图。我们可以看到在\(n_1\)固定的情况下，随着\(n_2\)的减小，F分布的图形也是变的越来越矮，相对也胖了一些，但是变化程度没有固定\(n_2\)的变化程度那么大。F分布并不是一个对称的图形，右侧有着一个长长的尾巴。

用Python绘画F分布概率密度函数图

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import f

if __name__ == '__main__':
    
    # 设置F分布的参数
    n1 = 20  # 自由度n1
    n2 = 40  # 自由度n2

    # 创建1000个0到50的等差数列
    x = np.linspace(0, 50, 1000)

    # 计算F分布的概率密度值
    f_pdf = f.pdf(x, n1, n2)

    # 绘制F分布的概率密度函数图
    plt.plot(x, f_pdf, color='blue', label='F({}, {})'.format(n1, n2))

    # 设置图表标题和轴标签
    plt.title('F Probability Density Function')
    plt.xlabel('x')
    plt.ylabel('f(x)')

    # 显示图例
    plt.legend()

    # 显示图表
    plt.show()

运行结果

正态总体的样本均值和样本方差的分布

设总体X的均值为μ，方差为\(σ^2\)，\(X_1,X_2,...,X_n\)是来自总体的一个样本，、\(S^2\)分别为样本均值和样本方差，则有

E()=μ，Var()=\(σ^2\over n\)，E(\(S^2\))=\(σ^2\)

这一组式子告诉我们，和\(S^2\)都是基于样本计算出来的统计量，而μ和\(σ^2\)却是总体的参数，之前我们说我们不可能拿到总体所包含的所有个体上的数据，所以我们也不可能去计算出总体的均值与方差，如果我们想了解总体的均值与方差，只能对总体进行观测，从总体中抽取一些个体，对个体进行观测，并从所得到的数据中去提取信息，基于这些信息来对总体的参数做出估计。

上面的这组等式为我们对总体的参数和样本的均值、方差之间建立了联系，使我们能够基于样本的信息去推测出总体的参数。这里需要注意的是，在这里我们并没有限定总体X的分布。

对E()=μ的证明

E()=\(E({X_1+X_2+...+X_n\over n})\)

\(={1\over n}E(X_1)+{1\over n}E(X_2)+...+{1\over n}E(X_n)\)

\(={1\over n}(μ+μ+...+μ)={1\over n}nμ=μ\)

这里因为\(X_1,X_2,...,X_n\)与总体X满足同分布，所以\(E(X_1),E(X_2),...,E(X_n)\)都是μ

对 Var()=\(σ^2\over n\)的证明

Var()=\(Var({X_1+X_2+...+X_n\over n})\)

\(={1\over n^2}Var(X_1+X_2+...+X_n)\)

\(={1\over n^2}Var(X_1)+{1\over n^2}Var(X_2)+...+{1\over n^2}Var(X_n)\)

\(={1\over n^2}(σ^2+σ^2+...+σ^2)={1\over n^2}nσ^2={σ^2\over n}\)

这里因为\(X_1,X_2,...,X_n\)与总体X满足同分布，所以\(Var(X_1),Var(X_2),...,Var(X_n)\)都是\(σ^2\)

• 定理一：\(X_1,X_2,...,X_n\)是来自正态总体N(μ,\(σ^2\))的一个样本，是样本均值，则有~N(μ,\(σ^2\over n\))
• 中心极线定理(Central Limit Theory)：中心极限定理说的是随机变量X的标准化函数

\(X-E(X)\over σ\)~N(0,1)

即该标准化函数近似服从标准正态分布，无论这个随机变量X服从何种分布。中心极限定理说明，在适当条件下，相互独立的随机变量之和经适当标准化之后，其分布近似于正态分布；不要求随机变量本身服从正态分布。

• 样本均值的分布

1、样本来自正态分布总体

总体服从标准正态分布 N (0,1)，从中抽取1000个样本(进行1000次抽样)，每次的样本容量=40，计算这1000个样本的样本均值，最终得到1000个样本均值的频数直方图如下

我们可以看到这1000个样本均值的分布非常接近于正态分布，上图中的红线代表理论均值0，绿线代表样本分布的均值，这两条线非常的接近，即样本分布的均值非常接近于0，这恰恰说明了样本均值=总体均值(期望)

2、样本来自非正态分布总体

总体服从均匀分布U(0,1)

这里a=0，b=1，\(1\over b-a\)为其概率密度函数。从中抽取1000个样本(进行1000次抽样)，每次的样本容量=40，计算这1000个样本的样本均值，最终得到1000个样本均值的频数直方图如下

我们可以看到这1000个样本均值的分布也非常接近于正态分布，如果我们增加样本个数的话，这个分布就可以更好的拟合正态分布。上图中的红线依然代表理论均值，绿线代表样本分布的均值，它们几乎也是重合的。

• 定理二：\(X_1,X_2,...,X_n\)是来自正态总体N(μ,\(σ^2\))的一个样本，、\(S^2\)分别为样本均值和样本方差，且相互独立，则有

\({(n-1)S^2\over σ^2}\)~\(χ^2(n-1)\)

该式将样本的方差与总体的方差建立的联系。

• 定理三： \(X_1,X_2,...,X_n\)是来自正态总体N(μ,\(σ^2\))的一个样本，、\(S^2\)分别为样本均值和样本方差，则有

\({mean(X)-μ\over {S\over \sqrt{n}}}\)~t(n-1)

这里的mean(X)即为，为样本均值。

由定理一，我们知道 ~N(μ,\(σ^2\over n\))，将该分布平移，伸缩，有

\({mean(X)-μ\over {σ\over \sqrt{n}}}\)～N(0,1)

我们将该式与定理三比较会发现左边的不同，一个是样本的标准差，服从自由度为n-1的t分布，一个是总体的标准差，服从标准正态分布。随着t分布自由度n的增大，我们知道它会逐渐逼近标准正态分布。这意味着当我们不知道总体的方差的时候，可以使用样本的方差来代替总体的方差，从而可以慢慢去逼近标准正态分布从而得到总体的方差。

• 定理四：(两个正态总体的样本均值和样本方差)

\(X_1,X_2,...,X_n\)是来自正态总体N(\(μ_1,σ_1^2\))的一个样本，、\(S_1^2\)分别为样本均值和样本方差；\(Y_1,Y_2,...,Y_n\)是来自正态总体N(\(μ_2,σ_2^2\))的一个样本，、\(S_2^2\)分别为样本均值和样本方差，则有

\({{S_1^2\over S_2^2}\over {σ_1^2\over σ_2^2}}\)~F(\(n_1-1,n_2-1\))

当\(σ_1^2=σ_2^2=σ^2\)时，

\({(mean(X)-mean(Y))-(μ_1-μ_2)\over S_W\sqrt{{1\over n_1}+{1\over n_2}}}\)~t(\(n_1+n_2-2\))
其中

\(S_W=\sqrt{{(n_1-1)S_1^2+(n_2-1)S_2^2}\over n_1+n_2-2}\)

编程理解中心极限定理

这里我们以之前的均匀分布为例来进行编程

import random
import matplotlib.pyplot as plt

def mean(data):
    return sum(data) / len(data)

if __name__ == '__main__':

    num_of_samples = 1000  # 样本数量
    sample_sz = 40  # 样本容量
    data = []  # 1000样本的均匀分布均值数据
    for _ in range(num_of_samples):
        data.append(mean([random.uniform(0.0, 1.0) for _ in range(sample_sz)]))
    plt.hist(data, bins='auto', rwidth=0.8)  # 样本均值的直方图
    plt.axvline(x=mean(data), c='red')  # 用红线画出样本均值的均值
    plt.show()

运行结果

如果我们觉得该图并没有那么的符合正态分布，那么我们可以增加样本数量

import random
import matplotlib.pyplot as plt

def mean(data):
    return sum(data) / len(data)

if __name__ == '__main__':

    num_of_samples = 100000  # 样本数量
    sample_sz = 40  # 样本容量
    data = []  # 100000样本的均匀分布均值数据
    for _ in range(num_of_samples):
        data.append(mean([random.uniform(0.0, 1.0) for _ in range(sample_sz)]))
    plt.hist(data, bins='auto', rwidth=0.8)  # 样本均值的直方图
    plt.axvline(x=mean(data), c='red')  # 用红线画出样本均值的均值
    plt.show()

这里我们将样本数量增加到了10万，运行结果

随机抽样、随机分配

样本是在总体中随机抽取的个体组成的，随着我们抽取出来的个体不同，我们就可以得到不同的样本。虽然这些样本是有变化的，但是好的样本的标准却是一致的——那就是我们希望我们的样本是具有代表性的，基于样本的结果可以很好的代表总体的结果，也就是说我们希望在样本上得到的结果也适用于在总体中没有被抽取到的样本个体的，这个过程就叫做随机抽样。

• 观察研究(Observational Study)

我们发现在我们周围有一些人非常热爱运动，而另外一些人从来不运动。现在我们有一个想法，我们认为运动与否会影响人的健康水平。

此时我们会从身边热爱运动的人中收集他们的健康水平的数据，然后再从身边不运动的人那里收集他们的健康水平数据。然后再将这两组人的健康水平进行对比，看看是否存在差异，从而作出运动与健康水平是否有关系的判断。此时我们得到的结果只能说明运动与健康水平是否相关，相关不代表因果，我们不能通过这种观察的研究去说明运动可以决定健康水平这种很强的因果的结论，因为在运动与否与健康水平之间可能有其他的混淆变量是我们所不知道的。

为了理解混淆变量，我们可以举这个例子，经过观察发现，相比于在结婚的时候没有办婚礼的人，举办了婚礼的人结婚以后离婚率更低。如果我们认为这个观察结果就是因果关系的话，那么就相当于在说办不办婚礼决定了离不离婚。但其实这两者之间之所以有关系，可能是因为第三个因素混淆变量在起作用，那就是这两个人是否是真爱，可能真爱的人更愿意花钱花精力去办一场婚礼。那么因此其实是真爱决定了两个人的离婚率的高低，而不是办不办婚礼决定的。

• 实验(Experiment)

实验可以去帮我们验证因果关系。我们首先先从总体中抽取一个样本，然后要把这个样本中的个体进行随机分配，把一半的个体分配到运动组中，把一半的个体分配到不运动组中。

对于运动组，我们会有计划的安排他们进行各种运动；对于不运动组，我们会要求他们在实验期间不进行任何运动。在一段时间后，我们就可以在这两组中收集他们的健康水平的数据，基于这两组数据，我们去进行比较和分析，如果这个时候我们依然发现运动组和不运动组在健康水平上有着明显的差异的话，我们就可以说运动与否决定了健康水平这个暗示了因果的结论，至少是适用于当前这个样本的。那么它能否推广到总体中去就取决于这个样本是否能够代表总体。

由以上我们可以得出以下的结论

参数估计

点估计

设总体X的分布函数的形式已知，但它的一个或多个参数未知，借助于总体X的一个样本来估计总体未知参数的值的问题称为参数的点估计问题。

• 估计量、估计值

θ是待估参数，\(X_1,X_2,...,X_n\)是总体X的一个样本，\(x_1,x_2,...,x_n\)是一个样本值。

点估计的问题就是要构造一个适当的统计量(估计量)，用它的观察值作为未知参数的近似值(估计值)。

这是我们之前取的正态分布总体和均匀分布总体的1000个样本的样本均值的频数直方图，上图直观的告诉了我们——对于不同的样本值，估计值一般是不同的。换句话说，总体的参数没有变，随着样本的变化，得到的样本值也会变化，基于样本值得到的估计值也会变化。这也就反映了我们的估计值具有不确定性，它跟真值之间是有一定的误差的。

• 估计量的评选标准

对于同一个参数，使用不同的估计方法求出的估计量可能不相同。使用哪个好？

1. 无偏性：若估计量的数学期望存在，并且该期望等于总体参数，则称为无偏估计。
2. 对于某些样本值，由这一估计量得到的估计值相对于真值来说偏大或偏小，但是反复将这一估计量使用多次，就"平均"来说其偏差为零。
3. E(估计值)-真值称为系统误差；无偏估计的实际意义就是无系统误差。

之前我们说过，设总体X的均值为μ，方差为\(σ^2\)，\(X_1,X_2,...,X_n\)是来自总体的一个样本，、\(S^2\)分别为样本均值和样本方差，则有

E()=μ，E(\(S^2\))=\(σ^2\)

说明不论总体服从什么分布，样本均值是总体均值的无偏估计；样本方差是总体方差的无偏估计。

而 Var()=\(σ^2\over n\)为样本均值的方差，表示样本抽样分布的离散程度，叫做标准误(standard error of mean;SE)

• 样本方差：除以n vs 除以n-1

在描述统计里，我们有

\(σ^2={\sum_{i=1}^n(X_i-mean(X))^2\over n}\)

它表示总体的方差

在样本方差中，我们有

\(S^2={\sum_{i=1}^n(X_i-mean(X))^2\over n-1}\)

它表示为总体方差的无偏估计。如果使用除以n的式子，则如果是基于样本计算的，则与总体方差有系统偏差。

我们以样例来说明，总体服从标准正态分布N(0,1)；样本容量=40；样本个数=1000，从而得到1000个样本方差的均值，如果使用\({\sum_{i=1}^n(X_i-mean(X))^2\over n}\)来计算每一个样本的方差，再画出所有样本方差的分布，就有

上图中的红线为总体的方差1，绿线为样本方差的均值，我们会看到绿线与红线之间有一定的距离。

 如果使用\({\sum_{i=1}^n(X_i-mean(X))^2\over n-1}\)来计算每一个样本的方差，再画出所有样本方差的分布，就有

上图中的绿线就非常的接近红线，虽然并未完全重合，但这个差距会随着样本容量的增加和样本个数的增加而逐渐减小直至消失。

import matplotlib.pyplot as plt
import numpy as np

def samle(num_of_samples, sample_sz):
    data = []
    for _ in range(num_of_samples):
        # 随机产生一组(40个)标准正态分布的数据
        x = np.random.normal(loc=0, scale=1, size=sample_sz)
        data.append(np.var(x))  # 获取这组数据的方差(\n方式)
    return data

if __name__ == '__main__':

    data = samle(1000, 40)
    plt.hist(data, bins='auto', rwidth=0.8)  # 画出1000个方差数据的直方图
    plt.axvline(x=np.mean(data), c='black')  # 画出这组方差的均值
    plt.axvline(x=1, c='red')  # 画出标准正态分布方差的均值
    plt.show()

运行结果

import matplotlib.pyplot as plt
import numpy as np

def samle(num_of_samples, sample_sz):
    data = []
    for _ in range(num_of_samples):
        # 随机产生一组(40个)标准正态分布的数据
        x = np.random.normal(loc=0, scale=1, size=sample_sz)
        data.append(np.var(x, ddof=1))  # 获取这组数据的方差(\(n-1)方式)
    return data

if __name__ == '__main__':

    data = samle(1000, 40)
    plt.hist(data, bins='auto', rwidth=0.8)  # 画出1000个方差数据的直方图
    plt.axvline(x=np.mean(data), c='black')  # 画出这组方差的均值
    plt.axvline(x=1, c='red')  # 画出标准正态分布方差的均值
    plt.show()

运行结果

总体服从均匀分布(0,1)；样本容量=40；样本个数=1000， 如果使用\({\sum_{i=1}^n(X_i-mean(X))^2\over n}\)来计算每一个样本的方差，再画出所有样本方差的分布，就有

如果使用\({\sum_{i=1}^n(X_i-mean(X))^2\over n-1}\)来计算每一个样本的方差，再画出所有样本方差的分布，就有

import matplotlib.pyplot as plt
import numpy as np

def samle(num_of_samples, sample_sz):
    data = []
    for _ in range(num_of_samples):
        # 随机产生一组(40个)0到1均匀分布的数据
        x = np.random.uniform(0, 1, sample_sz)
        data.append(np.var(x))  # 获取这组数据的方差(\n方式)
    return data

if __name__ == '__main__':

    data = samle(1000, 40)
    plt.hist(data, bins='auto', rwidth=0.8)  # 画出1000个方差数据的直方图
    plt.axvline(x=np.mean(data), c='black')  # 画出这组方差的均值
    plt.axvline(x=1 / 12, c='red')  # 画出0到1均匀分布方差的均值
    plt.show()

运行结果

import matplotlib.pyplot as plt
import numpy as np

def samle(num_of_samples, sample_sz):
    data = []
    for _ in range(num_of_samples):
        # 随机产生一组(40个)0到1均匀分布的数据
        x = np.random.uniform(0, 1, sample_sz)
        data.append(np.var(x, ddof=1))  # 获取这组数据的方差(\(n-1)方式)
    return data

if __name__ == '__main__':

    data = samle(1000, 40)
    plt.hist(data, bins='auto', rwidth=0.8)  # 画出1000个方差数据的直方图
    plt.axvline(x=np.mean(data), c='black')  # 画出这组方差的均值
    plt.axvline(x=1 / 12, c='red')  # 画出0到1均匀分布方差的均值
    plt.show()

运行结果

• 估计量的评选标准

对于同一个参数，使用不同的估计方法求出的估计量可能不相同。使用哪个好？

1. 有效性：有两个无偏估计\(θ_1\)和\(θ_2\)，如果在样本容量n相同的情况下，\(θ_1\)比\(θ_2\)更密集在真值附近，就认为\(θ_1\)比\(θ_2\)更理想；
2. 由于方差是随机变量取值与其数学期望的偏离程度的测量，所以无偏估计以方差最小者为好；
   1. 
   2. 可以用上图中的两个箭靶来说明，箭靶的中心就是真值，上面绿色的点就是环数。无论左靶还是右靶，它们绿点的均值都在靶心，但是右靶比左靶偏离靶心要更加严重，我们可以把这种偏离理解为方差。对于左右两个靶的均值都是无偏估计，那么显然左靶更好，因为它的绿点偏离靶心的程度更小。
3. 无偏性和有效性都是在样本容量n固定的前提下提出的。
4. 相合性：我们希望随着样本容量的增大，一个估计量的值稳定于代估参数的真值。满足此条件的估计量为相合估计量。

• 样本均值满足相合性

总体服从标准正态分布N(0,1)；样本容量=20，70，120，...，10000

上图中横轴为样本容量，纵轴为样本均值与总体均值的差。我们可以看到当n比较小的时候，震动的幅度还是很大的，随着样本容量的增加，震动的幅度越来越小，并且是围绕着0进行震动的。

• 样本方差满足相合性

总体服从标准正态分布N(0,1)；样本容量=20，70，120，...，10000

上图中横轴为样本容量，纵轴为样本方差与总体方差的差。我们可以看到当n比较小的时候，震动的幅度也是很大的，随着样本容量的增加，震动的幅度越来越小，并且是围绕着0进行震动的。

import numpy as np
import matplotlib.pyplot as plt

if __name__ == '__main__':

    sample_means = []
    sample_vars = []
    indices = []

    for sz in range(20, 10001, 50):
        indices.append(sz)
        # 随机产生sz个标准正态分布的数据
        sample = np.random.normal(loc=0, scale=1, size=sz)
        # 获取该组数据的均值
        sample_means.append(np.mean(sample))
        # 获取该组数据的方差
        sample_vars.append(np.var(sample, ddof=1))

    plt.plot(indices, sample_means)  # 画出所有组数据的均值折线图
    plt.plot(indices, sample_vars)  # 画出所有组数据的方差折线图
    plt.show()

运行结果

区间估计(interval estimate)

1. 对于一个未知量，我们在测量或计算时，不仅要得到近似值，还要估计误差，即近似值的精确程度/所求真值所在范围。
2. 对于未知参数，我们不仅要得到近似值(点估计)，还希望估计出一个范围(区间)，并希望知道这个范围包含参数真值的可信程度。这种形式的估计称为区间估计，这样的区间称为置信区间。

• 置信区间

设总体的分布函数含有一个未知参数θ，θ∈Θ(Θ是θ可能的取值范围)；对于给定值α(0<α<1)，若由来自X的样本\(X_1,X_2,...,X_n\)确定的两个统计量=(\(X_1,X_2,...,X_n\))和=(\(X_1,X_2,...,X_n\))(<)对于任意θ∈Θ满足

P{(\(X_1,X_2,...,X_n\))<θ<(\(X_1,X_2,...,X_n\))}≥1-α

则称随机区间(,)是θ的置信水平为1-α的置信区间，其中称为置信下限，称为置信上限，1-α称为置信水平。

固定样本容量n，若反复抽样多次，每个样本值确定一个区间，每个这样的区间要么包含θ的真值，要么不包含θ的真值。

上图中的虚线表示总体未知参数θ的真值，每一个竖线上的点代表点估计，每个点估计和它所在的区间都是从一个样本上得到的，所有的区间都对应一个样本。在图上我们可以看到，代表真值的虚线与垂直的蓝线有的是有交集的，有的是没有交集的(只有一个)，对于有交集的蓝线的区间就是包含了θ的真值，而单独的那个没有交集的蓝线的区间就是不包含θ真值的区间。对于置信区间要么包含真值，要么不包含真值，而不能说它以多大的概率包含真值。

按大数定律，在这么多区间中，包含真值的约占100*(1-α)%，不包含真值的占100*α%。如果α=0.05，那么95%的区间包含真值，5%的区间不包含真值。

• 求未知参数θ的置信区间的步骤

1. 寻求一个样本\(X_1,X_2,...,X_n\)和一个统计量W(\(X_1,X_2,...,X_n\);θ)，使统计量W的分布不依赖于θ和其他未知参数，统计量W的构造，通常可以从θ的点估计着手。
2. 对于给定的置信水平1-α，定出两个常数a和b，使得P(a<W(\(X_1,X_2,...,X_n\);θ)<b)=1-α，若能从中得到θ的不等式(\(X_1,X_2,...,X_n\))<θ<(\(X_1,X_2,...,X_n\))则(,)为θ的一个置信水平为1-α的置信区间。

设总体X~N(μ,\(σ^2\))，\(σ^2\)已知，μ为未知，设\(X_1,X_2,...,X_n\)是来自X的样本，求μ的置信水平为1-α的置信区间。

第1步， W(\(X_1,X_2,...,X_n\);θ)，W的分布不依赖于θ和其他未知参数，这里我们要构建的W=\({mean(X)-μ\over {σ\over \sqrt{n}}}\)~N(0,1)，这里的mean(X)，即就是样本均值的点估计，μ为总体的均值相当于统计量中的θ，分母中的σ和n都是已知数，那么该统计量W就服从标准正态分布。

第2步， P(a<W(\(X_1,X_2,...,X_n\);θ)<b)=1-α，由于W服从标准正态分布，我们来看一看标准正态分布

上图中的\(z_α\)为标准正态分布的上α分位点，含义为\(z_α\)往上部分的面积=α，用公式来表示就是\(P\{X>z_α\}=α\),(0<α<1)，这里的X可以替换成我们的W；\(-z_α\)为标准正态分布的下α分位点，含义为\(-z_α\)往下部分的面积=α，用公式来表示就是\(P\{X<-z_α\}=α\),(0<α<1)。那么从\(-z_α\)到\(z_α\)部分的面积就为1-2α。(这里的面积其实就是概率)

再回到P(a<W(\(X_1,X_2,...,X_n\);θ)<b)=1-α，我们就可以得出

在\(-z_{α\over 2}\)到\(z_{α\over 2}\)之间的面积就是1-α。自然而然，我们就可以得出\(a=-z_{α\over 2}\)，\(b=z_{α\over 2}\)。由此我们可以得出

\(P\{|{mean(X)-μ\over {σ\over \sqrt{n}}}|<z_{α\over 2}\}=1-α\)

它就是上图中中间区域的面积。解不等式，有

\(P\{mean(X)-{σ\over \sqrt{n}}z_{α\over 2}<μ<mean(X)+{σ\over \sqrt{n}}z_{α\over 2}\}=1-α\)

由此可得μ的置信区间为

\((mean(X)-{σ\over \sqrt{n}}z_{α\over 2},mean(X)+{σ\over \sqrt{n}}z_{α\over 2})\)

其中\(mean(X)-{σ\over \sqrt{n}}z_{α\over 2}\)为置信下限，\(mean(X)+{σ\over \sqrt{n}}z_{α\over 2}\)为置信上限。\(σ\over \sqrt{n}\)为标准误，它刻画的是样本均值离散的程度；\(z_{α\over 2}\)是一个关键值，它跟样本无关，只跟α和标准正态分布有关；\({σ\over \sqrt{n}}z_{α\over 2}\)叫做误差范围(margin or error;ME)。将该置信区间合写为

\((mean(X)(+/-){σ\over \sqrt{n}}z_{α\over 2})\)

我们重新来看这个真值与置信区间的图与置信区间的表达式一起来看看，上图中的点估计是随着样本均值mean(X)即变化而变化的，而区间的宽度是不变的，因为\({σ\over \sqrt{n}}z_{α\over 2}\)是不变的。假设现在只有一个样本，对置信区间的解读为：现在得到的区间属于那些包含μ的区间的可信程度为100*(1-α)%，或"该区间包含真值"这一陈述的可信程度为100*(1-α)%。

当然对于 P(a<W(\(X_1,X_2,...,X_n\);θ)<b)=1-α，对a,b的取值并不是唯一的，比如上图中，除了 \(-z_{α\over 2}\)到\(z_{α\over 2}\)之间的面积为1-α外，也有可能两条红线之间的面积或者两条绿线之间的面积都为1-α。因此置信水平为1-α的置信区间并不是唯一的，我们之所以选择了  \(-z_{α\over 2}\)到\(z_{α\over 2}\) 这组对称区间，是因为对于正态分布这种单峰且对称的情况，对称的区间长度最短，意味着估计的精度最高。比如说我们来估计程序员的收入水平，其中一个置信水平为95%的置信区间是2W-4W，另外一个置信水平为95%的置信区间是2W3到2W8，这个区间的长度显然小于2W-4W，并且这个区间的精度更高，它相当于给了我们一个更加窄的范围。因此我们在找置信区间的时候，我们希望找到那个精度最高的区间。所以在上面的例子中，我们求解的是这种对称区间。

• 置信区间与样本容量的关系

\((mean(X)(+/-){σ\over \sqrt{n}}z_{α\over 2})\)

1. 标准误\(σ\over \sqrt{n}\)会随着样本容量的增加而减小。
2. 误差范围\({σ\over \sqrt{n}}z_{α\over 2}\)会随着样本容量的增加而减小。

上图中的点依然是点估计，它们的数值都是一样的。我们可以看到随着n的增大，误差的范围越来越小，所以置信区间的宽度也就越来越小，也意味着精度越来越高。这也就告诉我们，样本容量是一个非常重要的指标，在我们进行数据分析的时候，希望使用的是大样本，在大样本上得到的估计(不管是点估计还是区间估计)的精度会更高，基于大样本得到的结论也会更加的可靠。

置信区间：一个正态总体的情况

• 方差已知，求均值的置信区间

示例1：歌曲的时长服从正态分布 X~N(μ,\(σ^2\))，\(σ^2\)=1；手机里有100首歌曲，这100首歌曲的平均时长是4分钟，求μ的置信水平为95%的置信区间。

由题意，我们可以知道

1. 样本均值=4
2. 样本容量n=100
3. 总体方差\(σ^2\)=1，标准差σ=1
4. α=0.05

由我们之前得到的置信区间有

\((mean(X)(+/-){σ\over \sqrt{n}}z_{α\over 2})\)

在上式中，我们唯一不知道的就是\(z_{α\over 2}\)，要得到该值，有两种方法，一种就是查表，另一种就是使用Python函数来得到，这里我们先看查表。

因为我们求的是置信水平为95%的置信区间，即中间的面积为0.95，那么两边的面积就为0.05，再两边各分一半，就是0.025，那么对应于均值的位置到\(z_{α\over 2}\)的面积就是0.475。我们现在要查的就是正态分布的Z分数表。

我们需要在上表中找到0.475这个数字，要得到该数值的Z分数，只要将该数值对应的行数与列数相加即可，则95%置信水平对应的Z值为1.9+0.06=1.96。

\((mean(X)(+/-){σ\over \sqrt{n}}z_{α\over 2})=(4(+/-){1\over \sqrt{100}}*1.96)=(3.804,4.196)\)

这就是我们得到的总体均值μ的置信水平为95%的置信区间。这一区间属于那些包含真值的区间的可信程度为95%，又可以说这一区间包含真值这一陈述的可信程度为95%，我们不能说该区间以95%的概率包含真值。

编程实现，我们的样本是18岁和35岁的人们薪资，求总体均值的置信区间

import numpy as np
from scipy.stats import norm
import math

def mean_ci_est(data, alpha, sigma):
    '''
    求均值的置信区间
    :param data: 样本
    :param alpha: 置信水平
    :param sigma: 总体标准差
    :return: 均值的置信区间
    '''
    n = len(data)  # 样本容量
    sample_mean = np.mean(data)  # 样本均值
    se = sigma / math.sqrt(n)  # 标准误
    z_value = abs(norm.ppf(alpha / 2))  # 正态分布的Z分数
    return sample_mean - se * z_value, sample_mean + se * z_value


if __name__ == '__main__':

    salary_18 = np.array([1484, 785, 1598, 1366, 1716, 1020, 1716, 785, 3113, 1601])
    salary_35 = np.array([902, 4508, 3809, 3923, 4276, 2065, 1601, 553, 3345, 2182])
    print(mean_ci_est(salary_18, 0.05, 1))
    print(mean_ci_est(salary_35, 0.05, 1))

运行结果

(1517.7802049676955, 1519.0197950323047)
(2715.7802049676957, 2717.0197950323045)

• 方差未知，求均值的置信区间

设总体X~N(μ,\(σ^2\))，\(σ^2\)未知，μ为未知，设\(X_1,X_2,...,X_n\)是来自X的样本，为样本均值，\(S^2\)为样本方差，求μ的置信水平为1-α的置信区间。

由于\(σ^2\)未知，所以我们不能使用 \({mean(X)-μ\over {σ\over \sqrt{n}}}\)~N(0,1)这个统计量，只能够使用样本方差来代替总体方差，则我们可以使用

W=\({mean(X)-μ\over {S\over \sqrt{n}}}\)~t(n-1)

这个统计量

 P(a<W(\(X_1,X_2,...,X_n\);θ)<b)=1-α

由上图可知，从\(-t_{α\over 2}(n-1)\)到\(t_{α\over 2}(n-1)\)之间的面积为1-α，则\(a=-t_{α\over 2}(n-1)\)，\(b=t_{α\over 2}(n-1)\)。由此我们可以得出

\(P\{|{mean(X)-μ\over {S\over \sqrt{n}}}|<t_{α\over 2}(n-1)\}=1-α\)

最终可得μ的置信区间为

\((mean(X)(+/-){S\over \sqrt{n}}t_{α\over 2}(n-1))\)

示例2： 歌曲的时长服从正态分布 X~N(μ,\(σ^2\))，\(σ^2\)未知；手机里有100首歌曲，这100首歌曲的平均时长是4分钟，方差为1.44，求μ的置信水平为95%的置信区间。

由题意，我们可以知道

1. 样本均值=4
2. 样本容量n=100
3. 样本方差\(S^2\)=1.44，标准差S=1.2
4. α=0.05

这里的置信区间为

\((mean(X)(+/-){S\over \sqrt{n}}t_{α\over 2}(n-1))\)

在上式中，我们唯一不知道的就是\(t_{α\over 2}(n-1)\)，这里我们直接给出该值，95%置信水平，自由度=99，对应的t值为1.98

\((mean(X)(+/-){S\over \sqrt{n}}t_{α\over 2}(n-1))=4(+/-){1.2\over \sqrt{100}}*1.98=(3.762,4.238)\)

这一区间属于那些包含真值的区间的可信程度为 95%，又可以说这一区间包含真值这一陈述的可信程度为 95%。

编程实现，我们的样本是18岁和35岁的人们薪资，求总体均值的置信区间

import numpy as np
from scipy.stats import t
import math

def mean_ci_est(data, alpha):
    '''
    求均值的置信区间
    :param data: 样本
    :param alpha: 置信水平
    :return: 均值的置信区间
    '''
    n = len(data)  # 样本容量
    sample_mean = np.mean(data)  # 样本均值
    s = np.std(data, ddof=1)  # 样本标准差
    se = s / math.sqrt(n)  # 标准误
    t_value = abs(t.ppf(alpha / 2, n - 1))  # t值
    return sample_mean - se * t_value, sample_mean + se * t_value


if __name__ == '__main__':

    salary_18 = np.array([1484, 785, 1598, 1366, 1716, 1020, 1716, 785, 3113, 1601])
    salary_35 = np.array([902, 4508, 3809, 3923, 4276, 2065, 1601, 553, 3345, 2182])
    print(mean_ci_est(salary_18, 0.05))
    print(mean_ci_est(salary_35, 0.05))

运行结果

(1042.5360403818845, 1994.2639596181157)
(1687.640703888853, 3745.1592961111473)

• 均值未知，求方差的置信区间

设总体X~N(μ,\(σ^2\))，\(σ^2\)未知，μ为未知，设\(X_1,X_2,...,X_n\)是来自X的样本，为样本均值，\(S^2\)为样本方差，求\(σ^2\)的置信水平为1-α的置信区间。

由于μ未知， 所以我们不能使用 \({mean(X)-μ\over {σ\over \sqrt{n}}}\)~N(0,1)这个统计量，也无法使用\({mean(X)-μ\over {S\over \sqrt{n}}}\)~t(n-1)这个统计量，我们可以使用

W=\({(n-1)S^2\over σ^2}\)~\(χ^2(n-1)\)

这个统计量

 P(a<W(\(X_1,X_2,...,X_n\);θ)<b)=1-α

虽然卡方分布的概率密度函数不是对称的，但我们仍在习惯上取对称的分位点，这里所谓的对称并不是位置上的对称，而是指它们分割的面积相等，这里都是\(α\over 2\)。由上图可知，从\(χ_{1-{α\over 2}}^2(n-1)\)到\(χ_{α\over 2}^2(n-1)\)之间的面积为1-α，则\(a=χ_{1-{α\over 2}}^2(n-1)\)，\(b=χ_{α\over 2}^2(n-1)\)。由此我们可以得出

\(P\{χ_{1-{α\over 2}}^2(n-1)<{(n-1)S^2\over σ^2}<χ_{α\over 2}^2(n-1)\}=1-α\)

最终可得\(σ^2\)的置信区间为

\(({(n-1)S^2\over χ_{α\over 2}^2(n-1)},{(n-1)S^2\over χ_{1-{α\over 2}}^2(n-1)})\)

示例3：歌曲的时长服从正态分布 X~N(μ,\(σ^2\))，μ未知；手机里有100首歌曲，这100首歌曲的时长方差为1.44，求\(σ^2\)的置信水平为95%的置信区间。

由题意，我们可以知道

1. 样本容量n=100
2. 样本方差\(S^2\)=1.44
3. α=0.05

这里的置信区间为

\(({(n-1)S^2\over χ_{α\over 2}^2(n-1)},{(n-1)S^2\over χ_{1-{α\over 2}}^2(n-1)})\)

在上式中，我们不知道的是\(χ_{α\over 2}^2(n-1)\)和\(χ_{1-{α\over 2}}^2(n-1)\)这两个值，我们直接给出该值，\(χ_{α\over 2}^2(n-1)=128.42\)，\(χ_{1-{α\over 2}}^2(n-1)=73.36\)

\(({(n-1)S^2\over χ_{α\over 2}^2(n-1)},{(n-1)S^2\over χ_{1-{α\over 2}}^2(n-1)})=({99*1.44\over 128.42},{99*1.44\over 73.36})=(1.11,1.94)\)

编程实现，我们的样本是18岁和35岁的人们薪资，求总体方差的置信区间

import numpy as np
from scipy.stats import chi2

def var_ci_est(data, alpha):
    '''
    求方差的置信区间
    :param data: 样本
    :param alpha: 置信水平
    :return: 方差的置信区间
    '''
    n = len(data)  # 样本容量
    s2 = np.var(data, ddof=1)  # 样本方差
    chi2_lower_value = chi2.ppf(alpha / 2, n - 1)  # 卡方值的上位点
    chi2_upper_value = chi2.ppf(1 - alpha / 2, n - 1)  # 卡方值的下位点
    return (n - 1) * s2 / chi2_upper_value, (n - 1) * s2 / chi2_lower_value


if __name__ == '__main__':

    salary_18 = np.array([1484, 785, 1598, 1366, 1716, 1020, 1716, 785, 3113, 1601])
    salary_35 = np.array([902, 4508, 3809, 3923, 4276, 2065, 1601, 553, 3345, 2182])
    print(var_ci_est(salary_18, 0.05))
    print(var_ci_est(salary_35, 0.05))

运行结果

(209357.6729819719, 1474810.356072325)
(978477.4012396415, 6892838.385031265)