动手学数理统计(1)

---

前言

前面的概率论专题，我们已经详细地学习了概率论的相关知识，并且针对某些特定问题使用概率论进行建模。但是在概率论的建模中，我们往往假设随机变量的分布已知，但是这个假设在实际中时很难做到的。现在，我们使用投资组合风险分析来说明：在前面的案例中，我们假设A、B股票的收益率服从$N(0.1,0.01)$和$N(0.3,0.04)$，但是现实时市场上并不会告诉我们股票的收益分布是什么，我们能做的只有根据历史数据进行合理的猜测与估算，而不同的股票的收益率分布都不一样。那么统计学能做的就是对不同的股票合理估算和猜测相应的收益率的分布，找到最佳的参数，如正态分布的$\mu$和${\sigma }^2$。不仅如此，统计学关心的不只有$\mu$和${\sigma }^2$是多少，也会关心比如$\mu$和${\sigma }^2$是否符合符合相关要求等等。

使用一个更通俗的例子来说，用于概率论的建模，现在假设中国人的身高近似服从正态分布，那究竟这个正态分布中的$\mu$和${\sigma }^2$的值是多少就是统计学做的事情，而这个工作也叫做参数估计，后面的内容会详细介绍。

一般认为，统计学是收集、分析、表述和解释数据的科学，统计学是一门处理数据的方法和技术的学科。

---

一、总体和样本

在一个统计问题中, 我们把研究对象的全体称为总体， 构成总体的每个成员称为个体。对大多数实际问题, 总体中的个体是一些实在的人或物，如：假如我们想要研究中国人的身高情况，那么全体中国人构成该研究问题的总体，而每一个中国人就是构成该研究问题的个体。但是，在实际的操作中，我们只关心每个中国人的身高，并不关心每个中国人的年龄、性别等。因此，每个中国人的数量指标————身高就变成了个体，而全体中国人的所有身高也就构成了总体。因此总体和个体就是一堆数字，这堆数字有大有小，每个数字出现的可能性又不太一样，看似有某个规律。因此，我们使用一个概率分布描述和归纳总体是合适的且合理的，换句话来说：总体就是一个概率分布，总体的数量指标就是服从该概率分布的一个随机变量。因此，从总体中抽样与从某个分布中抽样的含义是一样的。

举个例子：去银行排队办理业务，队伍中排了非常多的人，我们想要研究队伍中每个人的排队等待的时间情况。因此，这个问题的研究总体为队伍里排队的所有人，个体是每一个排队等待的人，我们想要关心的数量指标是排队的时间，有人排5min，有人排3.2min等等，因此队伍里排队的所有人的一个数量指标：等待时间构成了一个随机变量，假设这个随机变量服从参数为$\lambda$的指数分布$Exp(\lambda )$。但是分布中的参数$\lambda$是未知的，$\lambda$的大小决定了银行服务系统的质量，它直接影响银行的声誉。

本例中总体分布的类型是明确的, 是指数分布, 但是总体含有末知参数 $\lambda$, 因此总体不是一个特定的泊松分布。 要确定最终的总体分布, 就是要确定符合具体某个队伍的 $\lambda$， 这是统计学科的任务之一。

一般来说，总体分为：有限总体和无限总体, 大多数我们说的总体是无限总体。

由于总体是无限的，又或者说总体的个体数量太多，我们如果对每一个个体的数量指标都进行研究所花费的成本将是十分巨大的。为了了解总体的分布， 我们从总体中随机地抽取 $n$ 个个体， 记其指标值为 $x_1,x_2,\cdots$, $x_n$, 则 $x_1,x_2,\cdots ,x_n$ 称为总体的一个样本， $n$ 称为样本容量， 或简称样本量， 样本中的个体称为样品。

下面这段描述样本的概念是非常重要的，希望大家牢牢谨记！！！

样本具有所谓的二重性：一方面, 由于样本是从总体中随机抽取的， 抽取前无法预知它们的数值， 因此， 样本是一个随机变量， 用大写字母 $X_1,X_2,\cdots ,X_n$ 表示；另一方面， 样本在抽取以后经观测就有确定的观测值， 因此， 样本又是一组数值， 此时用小写字母 $x_1,x_2,\cdots ,x_n$ 表示是合适的。 为了描述的简单，我们只用小写字母表示样本$x_1,x_2,\cdots ,x_n$，不管样本是随机变量还是具体的数值，都用小写字母表示。

举个例子：某个肉类产品的规定净重为1kg，由于生产过程中具有随机性，不可能让所有实际生产的肉类品都是1kg，现在从同一批次生产的肉类产品随机抽取10份肉类品进行质量检查并称重，称重结果如下：
$$\begin{array}{llllllllll}0.98& 0.99& 0.95& 1.01& 1.04& 0.96& 0.96& 1.03& 0.98& 1.00\end{array}$$
这是一个容量为 10 的样本的观测值，对应的总体为该厂生产的肉类产品的净重。

从总体中抽取样本的方式多种多样，随心所欲，如：从大到小抽取、随机抽取等等。 为了能由样本对总体作出较可靠的推断, 就希望样本能很好地代表总体. 需要对抽样的方法提出一些要求,最常用的抽样方法就是“简单随机抽样”，简单随机抽样需要满足两个要求：

• 从总体中抽取的样本具有代表性：具有代表性要求总体中每一个个体都有同等机会被选入样本中，也就意味着样本中的每一个样品$x_i$与总体X有相同的分布，简称“同分布”。
• 从总体中抽取的样本具有独立性：具有独立性即要求样本中每一样品的取值不影响其他样品的取值， 也就意味着 $x_1,x_2,\cdots ,x_n$ 之间相互独立。

总结起来：在简单随机抽样这种抽样方法下，样本中的每一个样品$x_1,x_2,\cdots ,x_n$之间独立同分布，同分布于总体分布，简称：iid。用简单随机抽样方法得到的样本称为简单随机样本，也简称样本。

设总体 $X$ 具有分布函数 $F(x),x_1,x_2,\cdots ,x_n$ 为取自该总体的容量为 $n$ 的样本，则样本联合分布函数为
$$F\left(x_1,x_2,\cdots ,x_n\right)=\prod _{i=1}^{n}F\left(x_i\right).$$
这个公式在极大似然估计的时候会再次碰到，是一条十分重要的公式！

有没有情况是不符合独立同分布的假设呢？其实很多情况下是不满足独立同分布的假设的，比如：时间序列数据，前后之间存在关联，因此不是相互独立的。

二、经验分布函数与直方图

（1）经验分布函数：

统计学的一个重要核心就是使用样本信息估计总体信息，有的时候总体往往是未知的，我们只能通过多次试验的样本（即实际值）来推断总体。经验分布函数就是使用样本信息构造的分布函数近似未知的总体分布函数：

设 $x_1,x_2,\cdots ,x_n$ 是取自总体分布函数为 $F(x)$ 的样本， 若将样本观测值由小到大进行排列， 记为 $x_{(1)},x_{(2)},\cdots ,x_{(n)}$, 则 $x_{(1)},x_{(2)},\cdots ,x_{(n)}$ 称为有序样本， 用有序样本定义如下函数
$$F_n(x)=\{\begin{array}{ll}0,& \text{ 当 }x<x_{(1)},\\ k/n,& \text{ 当 }{x}_{(k)}⩽x<x_{(k+1)},k=1,2,\cdots ,n-1,\\ 1,& \text{ 当 }x⩾x_{(n)},\end{array}$$
则 $F_n(x)$ 是一非减右连续函数， 且满足
$$F_n(-\infty )=0\text{ 和 }{F}_n(\infty )=1.$$
由此可见， $F_n(x)$ 是一个分布函数， 称 $F_n(x)$ 为该样本的经验分布函数。

为什么样本信息可以用来推断总体呢？格纹利柯定理给我们指明了一个方向：

设 $x_1,x_2,\cdots ,x_n$ 是取自总体分布函数为 $F(x)$ 的样本， $F_n(x)$ 是其经验分布函数， 当 $n\to \infty$ 时， 有
$$P\left(\underset{-\infty <x<\infty }{\sup }\left|F_n(x)-F(x)\right|\to 0\right)=1$$
当 $n$ 相当大时，经验分布函数是总体分布函数 $F(x)$ 的一个良好的近似。 因此，在经典统计学中一切统计推断都是以样本为依据推断总体。

例子：随机观察总体X, 得到一个容量为 10 的样本:
$$3.2,2.5,-2,2.5,0,3,2,2.5,2,4$$
求 $\mathrm{X}$ 经验分布函数。

解：首先，将样本排序 $-2<0<2=2<2.5=2.5=2.5<3<3.2<4$

于是，根据公式可以得到经验分布函数：
$$F_{10}(x)=\{\begin{array}{cc}0,& x<-2\\ 1/10,& -2\le x<0\\ 2/10,& 0\le x<2\\ 4/10,& 2\le x<2.5\\ 7/10,& 2.5\le x<3\\ 8/10,& 3\le x<3.2\\ 9/10,& 3.2\le x<4\\ 1,& x\ge 4\end{array},$$

（2）直方图：频数直方图和频率直方图

直方图是数值数据分布的精确图形表示， 这是一个连续变量（定量变量）的概率分布的估计，被卡尔·皮尔逊（Karl Pearson）首先引入。

# 引入相关工具库
import numpy as np 
import pandas as pd 
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use("ggplot")
import warnings 
warnings.filterwarnings("ignore")
plt.rcParams['font.sans-serif']=['SimHei','Songti SC','STFangsong']
plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
import seaborn as sns 

# 频数直方图
x_samples = np.random.randn(1000)
plt.hist(x_samples, bins=10,color='green',alpha=0.6)  # bins=10代表10根柱子
plt.xlabel("x")
plt.ylabel("频数 n")
plt.title("频数直方图")
plt.show()

# 频率直方图
x_samples = np.random.randn(1000)
plt.hist(x_samples, bins=10,color='blue',alpha=0.6,density=True)  # bins=10代表10根柱子
plt.xlabel("x")
plt.ylabel("频率 p")
plt.title("频率直方图")
plt.show()

三、统计量与三大抽样分布

（1）统计量：

在统计学中，我们希望使用样本信息推断总体，而且样本也来源于总体，但是样本的信息较为分散，有时候还会显得杂乱无章。为了将这些信息分散的样本中有关总体的那部分信息汇集起来，从而反映总体的某些特征，需要对样本进行加工。前面的直方图就是一个很好的加工方式，但是直方图这样的图形化方式有一个致命的缺点就是无法数值化进行研究（例如比较两个总体某个特征的大小）。进行数值化最有效的方式就是构造关于样本的函数，函数是数值化最直接的办法，而且不同的关于样本的函数反映了总体的不同特征。统计量就是在这样的背景中诞生的：

设 $x_1,x_2,\cdots ,x_n$ 为取自某总体的样本， 若样本函数 $T=T\left(x_1,x_2,\cdots ,x_n\right)$ 中不含有任何末知参数， 则称 $T$ 为统计量。统计量的分布称为抽样分布。

例如：若 $x_1,x_2,\cdots ,x_n$ 为样本， 则 ${\sum }_{i=1}^n{x}_i,{\sum }_{i=1}^n{x}_i^2$ 以及 $F_n(x)$ 都是统计量，我们中学开始学习的样本平均数$\bar{x}=\frac{x_1+x_2+\cdots +x_n}{n}=\frac{1}{n}{\sum }_{i=1}^n{x}_i$和样本方差$s_n^2=\frac{1}{n}{\sum }_{i=1}^n{\left(x_i-\bar{x}\right)}^2$都是统计量，因为这些函数都是只关于样本的函数，不含未知参数。而当 $\mu ,{\sigma }^2$ 未知时， $x_1-\mu ,x_1/\sigma$ 等均不是统计量。值得注意的是：统计量由样本决定，从而统计量因样本而异，对于同一总体，由于抽取样本是具有随机性的，因此抽取不同的样本，统计量就不同，从而统计量也是一个随机变量。统计量的分布称为抽样分布。虽然统计量不依赖于任何参数，但统计量的分布一般依赖于未知参数。

下面介绍几个常用的统计量及其对应的抽样分布：

• （1.1）样本均值：

设 $x_1,x_2,\cdots ,x_n$ 为取自某总体的样本， 其算术平均值称为样本均值，一 般用 $\bar{x}$ 表示，即
$$\bar{x}=\frac{x_1+x_2+\cdots +x_n}{n}=\frac{1}{n}\sum _{i=1}^n{x}_i$$
如果把样本中的数据与样本均值的差称为偏差， 则样本所有偏差之和为 0， 即 ${\sum }_{i=1}^n\left(x_i-\bar{x}\right)=0$。

# 从总体/总体的分布中抽取样本并计算样本均值和计算偏差
## （1）从总体中抽取样本
X = np.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20])   # 假设总体为X
x_i = np.random.choice(X, 10, replace=False)  # 从总体X中抽取10个样本
x_mean = np.mean(x_i) # 计算样本均值
x_bias = np.sum(x_i-x_mean)  # 计算偏差和
print("样本均值为：",x_mean)
print("偏差和为：",x_bias)

样本均值为： 10.1
偏差和为： 3.552713678800501e-15

## （2）从总体分布中抽取样本，假设总体分布为N（0，1）
x_i = np.random.randn(10)  # 从总体分布N(0,1)中抽取10个样本
x_mean = np.mean(x_i) # 计算样本均值
x_bias = np.sum(x_i-x_mean)  # 计算偏差和
print("样本均值为：",x_mean)
print("偏差和为：",x_bias)

样本均值为： 0.06436695925574631
偏差和为： 4.440892098500626e-16

接下来，我们来看看样本均值的分布是什么？我们用一个实验给大家说明：

• （a）假设现在有20个数组成的总体；
• （b）每次从总体中抽取5个样本，计算样本均值；
• （c）重复b步骤5次，10次，20次，100次，1000次，10000次，100000次，观察不同重复次数的抽样的样本均值的分布情况。

def GetSampleDist(n, X):
    x_mean_list = []
    for i in range(n):
        x_i = np.random.choice(X, 5, replace=False)
        x_mean = np.mean(x_i)
        x_mean_list.append(x_mean)
    plt.hist(x_mean_list,color='orange',alpha=0.6,density=True)
    plt.xlabel("x")
    plt.ylabel("频率 p")
    plt.title("n="+str(n))
    plt.show()

X = np.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20])   # 假设总体为X, size=20
GetSampleDist(5, X)
GetSampleDist(10, X)
GetSampleDist(20, X)
GetSampleDist(100, X)
GetSampleDist(1000, X)
GetSampleDist(10000, X)
GetSampleDist(100000, X)

因此，从上面的实验可以看到，统计量—样本均值的分布，即样本均值的抽样分布当n越来越大时近似服从正态分布，具体来说：

设 $x_1,x_2,\cdots ,x_n$ 是来自某个总体的样本， $\bar{x}$ 为样本均值。

（1）若总体分布为 $N\left(\mu ,{\sigma }^2\right)$， 则 $\bar{x}$ 的精确分布为 $N\left(\mu ,{\sigma }^2/n\right)$;

（2） 若总体分布末知或不是正态分布， $E(x)=\mu ,\mathrm{Var}(x)={\sigma }^2$ 存在， 则 $n$ 较大时 $\bar{x}$ 的渐近分布为 $N\left(\mu ,{\sigma }^2/n\right)$。 这里渐近分布是指 $n$ 较大时的近似分布。

• （1.2）样本方差与样本标准差：

设 $x_1,x_2,\cdots ,x_n$ 为取自某总体的样本，则它关于样本均值 $\bar{x}$ 的平均偏差平方和
$$s_n^2=\frac{1}{n}\sum _{i=1}^n{\left(x_i-\bar{x}\right)}^2$$
称为样本方差。

样本标准差就是样本方差的算术平方根，即：$s_n=\sqrt{s_n^2}$。

在实际的应用中，我们往往不会直接使用样本方差，更倾向于使用样本标准差，因为样本方差的量纲（单位）与样本均值不一致，无法与样本均值进行加减运算。实际上，样本标准差可以与均值发生计算关系，表示数据的范围，如：$(\bar{x}-3s_n,\bar{x}+3s_n)$表示数据的范围在样本均值的三个标准差范围，而方差由于量纲不一致做不到与样本均值的联动。

另外，样本方差除了上述的表达式，还有另一个表达式：$s^2=\frac{1}{n-1}{\sum }_{i=1}^n{\left(x_i-\bar{x}\right)}^2$。$s^2$常被称为无偏方差，什么是无偏将在后面的内容讲解，这里大家只需要知道的就是：当样本量较大时，$s_n^2$与$s^2$相差不大，可以随意使用，当样本量较小时，计算样本方差最好使用无偏样本方差$s^2$。

值得注意的是：后面我们所说的样本方差都是指无偏样本方差$s^2$而不是$s_n^2$，请大家注意。

# 从总体/总体的分布中抽取样本并计算样本方差与样本标准差
## （1）从总体中抽取样本
X = np.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20])   # 假设总体为X
x_i = np.random.choice(X, 10, replace=False)  # 从总体X中抽取10个样本
x_sn2 = np.var(x_i,ddof=0) #样本方差
x_s2 = np.var(x_i,ddof=1) # 无偏样本方差
x_sn = np.std(x_i,ddof=0) # 样本标准差
x_s = np.std(x_i,ddof=1) # 无偏样本标准差
print("样本方差sn^2为：",x_sn2)
print("样本方差s^2为：",x_s2)
print("样本标准差sn^2为：",x_sn)
print("样本标准差s^2为：",x_s)

样本方差sn^2为： 38.64
样本方差s^2为： 42.93333333333333
样本标准差sn^2为： 6.216108107168021
样本标准差s^2为： 6.552353266829661

## （2）从总体分布中抽取样本，假设总体分布为N（0，1）
x_i = np.random.randn(10)  # 从总体分布N(0,1)中抽取10个样本
x_sn2 = np.var(x_i,ddof=0) #样本方差
x_s2 = np.var(x_i,ddof=1) # 无偏样本方差
x_sn = np.std(x_i,ddof=0) # 样本标准差
x_s = np.std(x_i,ddof=1) # 无偏样本标准差
print("样本方差sn^2为：",x_sn2)
print("样本方差s^2为：",x_s2)
print("样本标准差sn^2为：",x_sn)
print("样本标准差s^2为：",x_s)

样本方差sn^2为： 1.6432573626825586
样本方差s^2为： 1.825841514091732
样本标准差sn^2为： 1.2818960030683295
样本标准差s^2为： 1.3512370310540383

样本方差/样本标准差的分布：

事实上，样本方差/样本标准差的分布并没有像样本均值那样有完美的渐近分布，不同总体下的样本方差/样本标准差的分布都不一样，我们在下一节内容给大家介绍当总体服从正态分布条件下的与样本方差/样本标准差有关的分布，他是三大抽样分布中著名的卡方分布。在这里，我们可以发现如下规律：

设总体为 $X$ 方差存在， 即 $E(X)=\mu ,\mathrm{Var}(X)={\sigma }^2<\infty ,x_1,x_2,\cdots ,x_n$ 为 从该总体得到的样本， $\bar{x}$ 和 $s^2$ 分别是样本均值和样本方差, 则
$$\begin{array}{c}E(\bar{x})=\mu ,\mathrm{Var}(\bar{x})={\sigma }^2/n,\\ E\left(s^2\right)={\sigma }^2.\end{array}$$
此定理表明， 样本均值的期望与总体均值相同， 而样本均值的方差是总体方差的 $1/n$。

下面，我们使用一个实验来验证以上的观点：

• （a）假设总体为$X$，为了方便起见，$X～N（0，1）$；
• （b）对$X$抽取10个样本，计算10个样本的样本均值；
• （c）重复（b）步骤10000次，计算这10000个样本均值的样本均值和样本方差；

x_mean_list = []
for i in range(10000):
    x_i = np.random.randn(10)
    x_mean_list.append(np.mean(x_i))
print("标准正态分布的均值和方差为：",0,1)
print("1000个样本均值的样本均值为：",np.mean(x_mean_list))
print("1000个样本均值的样本方差为：",np.var(x_mean_list,ddof=1))

标准正态分布的均值和方差为： 0 1
1000个样本均值的样本均值为： 0.0004640811696833645
1000个样本均值的样本方差为： 0.09955649928109517

• （1.3）次序统计量及其分布（了解）

设 $x_1,x_2,\cdots ,x_n$ 是取自总体 $X$ 的样本， $x_{(i)}$ 称为该样本的第 $i$ 个次序统计量， 它的取值是将样本观测值由小到大排列后得到的第 $i$ 个观测值。 其中 $x_{(1)}=$ min ⁡ { x 1 , x 2 , ⋯   , x n } \min \left\{x_{1}, x_{2}, \cdots, x_{n}\right\} min{x1​,x2​,⋯,xn​} 称为该样本的最小次序统计量， x ( n ) = max ⁡ { x 1 , x 2 , ⋯   , x n } x_{(n)}=\max \left\{x_{1}, x_{2}, \cdots, x_{n}\right\} x(n)​=max{x1​,x2​,⋯,xn​} 称为该样本的最大次序统计量。

以上的定义十分难理解，我们使用一个实验说明什么是次序统计量，特别强调最小、最大次序统计量：

（a）最小次序统计量及其分布：

• （a.1）假设总体X服从以下分布：
  $$\begin{array}{cccc}x& 1& 2& 3\\ p& 1/3& 1/3& 1/3\end{array}$$
• （a.2）现在从总体中抽取10个样本，每次取10个样本中最小的那个样品值为最小次序统计量。
• （a.3）重复（a.2）步骤5次，10次，20次，100次，1000次，10000次，观察$n$越来越大时的最小次序统计量的稳定后分布。

理论计算的最小次序统计量的分布：
$$\begin{array}{cccc}{x}_{(1)}& 1& 2& 3\\ p& \frac{19}{27}=0.703& \frac{7}{27}=0.259& \frac{1}{27}=0.038\end{array}$$
（b）最小次序统计量及其分布：

• （b.1）假设总体$X$服从以下分布：
  $$\begin{array}{cccc}x& 1& 2& 3\\ p& 1/3& 1/3& 1/3\end{array}$$
• （b.2）现在从总体中抽取10个样本，每次取10个样本中最大的那个样品值为最大次序统计量。
• （b.3）重复（a.2）步骤5次，10次，20次，100次，1000次，10000次，观察n越来越大时的最大次序统计量的稳定后分布。

理论计算的最大次序统计量的分布：
$$\begin{array}{cccc}{x}_{(3)}& 1& 2& 3\\ p& \frac{1}{27}=0.038& \frac{7}{27}=0.259& \frac{19}{27}=0.703\end{array}$$

# 最小次序统计量及其分布
from scipy.stats import rv_discrete  # 自定义离散分布
x_k = np.arange(3)+1
p_k = np.array([1/3]*3)
X = rv_discrete(name='min', values=(x_k, p_k))
def Get_Min_Dist(n, X):
    Min_list = []
    for i in range(n):
        min_value = np.min(X.rvs(size=3))
        Min_list.append(min_value)
    # 统计每个值出现的次数
    xk_count = []
    for v in x_k:
        xk_count.append(np.sum(Min_list==v))
    # 画图
    plt.bar(x_k,np.array(xk_count)/n,color='blue',alpha=0.6)
    plt.xlabel("x")
    plt.ylabel("频率 p")
    plt.xlim(0,4)
    plt.title("n="+str(n))
    plt.show()

Get_Min_Dist(5, X)
Get_Min_Dist(10, X)
Get_Min_Dist(20, X)
Get_Min_Dist(100, X)
Get_Min_Dist(1000, X)
Get_Min_Dist(10000, X)

# 最大次序统计量及其分布
from scipy.stats import rv_discrete  # 自定义离散分布
x_k = np.arange(3)+1
p_k = np.array([1/3]*3)
X = rv_discrete(name='min', values=(x_k, p_k))
def Get_Max_Dist(n, X):
    Max_list = []
    for i in range(n):
        max_value = np.max(X.rvs(size=3))
        Max_list.append(max_value)
    # 统计每个值出现的次数
    xk_count = []
    for v in x_k:
        xk_count.append(np.sum(Max_list==v))
    # 画图
    plt.bar(x_k,np.array(xk_count)/n,color='blue',alpha=0.6)
    plt.xlabel("x")
    plt.ylabel("频率 p")
    plt.xlim(0,4)
    plt.title("n="+str(n))
    plt.show()

Get_Max_Dist(5, X)
Get_Max_Dist(10, X)
Get_Max_Dist(20, X)
Get_Max_Dist(100, X)
Get_Max_Dist(1000, X)
Get_Max_Dist(10000, X)

• （1.4）样本分位数与样本中位数及其抽样分布

在概率论中，我们知道中位数是指累计概率$p=0.5$时对应的随机变量值$x$，而样本抽样后的一组离散个数的取值，因此，我们定义样本中位数只需要查找抽样后的样本中排在中间的那个样品，具体来说：
$$m_{0.5}=\{\begin{array}{ll}x\left(\frac{n+1}{2}\right),& n\text{ 为奇数, }\\ \frac{1}{2}\left(\frac{x}{2}\right)+x\left(\frac{n}{2}+1\right)),& n\text{ 为偶数. }\end{array}$$
若 $n=5$， 则 $m_{0.5}=x_{(3)}$， 若 $n=6$, 则 $m_{0.5}=\frac{1}{2}\left(x_{(3)}+x_{(4)}\right)$。

# 从总体/总体的分布中抽取样本并计算样本样本中位数
## （1）从总体中抽取样本
X = np.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20])   # 假设总体为X
x_i = np.random.choice(X, 10, replace=False)  # 从总体X中抽取10个样本
x_mid = np.median(x_i) # 计算样本中位数
print("样本中位数为：",x_mid)

样本中位数为： 9.5

## （2）从总体分布中抽取样本，假设总体分布为N（0，1）
x_i = np.random.randn(10)  # 从总体分布N(0,1)中抽取10个样本
x_mid = np.median(x_i) # 计算样本中位数
print("样本中位数为：",x_mid)

样本中位数为： -0.12095156770978058

接下来，我们来看看样本p分位数是什么？

样本 $p$ 分位数 $m_p$ 可如下定义:
$$m_p=\{\begin{array}{ll}{x}_{([np+1])},& \text{ 若 }np\text{ 不是整数, }\\ \frac{1}{2}\left(x_{(np)}+x_{(np+1)}\right),& \text{ 若 }np\text{ 是整数. }\end{array}$$
用普通的语言理解就是：先将抽样的$n$个样本按从小到大排序，样本$p$分位数就是排在第$n\times p$位的样品值，当然有可能$n\times p$不是整数，按上述公式处理。

# 从总体/总体的分布中抽取样本并计算样本样本中位数
## （1）从总体中抽取样本
X = np.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20])   # 假设总体为X
x_i = np.random.choice(X, 10, replace=False)  # 从总体X中抽取10个样本
x_low = np.percentile(x_i,25) # 计算样本下四分位数
x_high = np.percentile(x_i,75) # 计算样本上四分位数
print("样本下四分位数为：",x_low)
print("样本上四分位数为：",x_high)

样本下四分位数为： 4.25
样本上四分位数为： 14.25

接下来，我们使用一组实验观察下样本p分位数/样本中位数的抽样分布：

• （a）假设总体为$X$，为了方便起见，$X$取1～20；
• （b）对$X$抽取10个样本，计算10个样本的样本中位数；
• （c）重复（b）步骤10次，20次，100次，1000次，10000次，100000次，观察样本中位数的分布；

def GetMidDist(n, X):
    x_mid_list = []
    for i in range(n):
        x_i = np.random.choice(X, 10, replace=False)
        x_mid = np.median(x_i)
        x_mid_list.append(x_mid)
    plt.hist(x_mid_list,color='blue',alpha=0.6,density=True)
    plt.xlabel("x")
    plt.ylabel("频率 p")
    plt.title("n="+str(n))
    plt.show()

X = np.array([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20])   # 假设总体为X, size=20
GetMidDist(5, X)
GetMidDist(10, X)
GetMidDist(20, X)
GetMidDist(100, X)
GetMidDist(1000, X)
GetMidDist(10000, X)
GetMidDist(100000, X)

因此，样本分位数的抽样分布的渐近分布为正态分布，当n越来越大时。具体来说：

设总体密度函数为 $p(x),x_p$ 为其 $p$ 分位数， $p(x)$ 在 $x_p$ 处连续且 $p\left(x_p\right)>$ 0 ， 则当 $n\to \infty$ 时样本 $p$ 分位数 $m_p$ 的渐近分布为
$$N\left(x_p,\frac{p(1-p)}{n\cdot p^2\left(x_p\right)}\right)$$
特别地, 对于样本中位数来说, 当 $n\to \infty$ 时有
$$N\left(x_{0.5},\frac{1}{4n\cdot p^2\left(x_{0.5}\right)}\right)$$

（2）三大抽样分布：

以上关于统计量及其分布的介绍，我们已经学习了几种常用的统计量及其分布，包括：样本均值、样本方差/样本标准差以及分位数/中位数。在这些统计量的构造以及抽样分布的探索上，我们并没有对总体分布的形式进行限定，这是一件好事，也是一件不好的事，因为如果对总体的分布做限定，会得出特定条件下的统计量及其抽样分布，但是会降低使用的广泛性。实际上，有很多实际问题都是以标准正态分布的前提下讨论的，因此基于此，讨论以标准正态分布为总体而构造的统计量及其抽样分布具有广泛的应用。以标准正态分布为总体构造的三个著名统计量：卡方统计量、F统计量与t统计量在实际中有广泛的用途。下面，我们来一一学习！

• （2.1）卡方统计量与卡方分布：

假设$x_1,x_2,...x_n$是标准正态分布$N(0,1)$为总体抽样的得到的样本（$x_1,x_2,...x_n$独立同分布于$N(0,1)$），则
$${\chi }^2=x_1^2+x_2^2+\cdots +x_n^2$$
的分布为自由度为n的${\chi }^2$分布，简称${\chi }^2\sim {\chi }^2(n)$，${\chi }^2$分布的密度函数为：
$$p(y)=\frac{(1/2{)}^{\frac{n}{2}}}{\Gamma (n/2)}{y}^{\frac{n}{2}-1}{\mathrm{e}}^{-\frac{y}{2}},y>0$$
可以用一句话记住卡方分布：$n$个标准正态分布的平方和服从自由度为$n$的卡方分布。

有同学可能会问，什么是自由度，这个概念不用深究，可以简单理解为可以自由变化的变量个数。有一个例子很好地阐述了自由度的概念：在无偏样本方差$s^2=\frac{1}{n-1}{\sum }_{i=1}^n{\left(x_i-\bar{x}\right)}^2$中，为什么分母是$n-1$而不是$n$呢？是因为在这个式子中只有$n-1$个可以自由变化的变量数，$n$个样本，其中在计算样本均值时需要花费一个方程，即$\bar{x}=\frac{x_1+x_2+\cdots +x_n}{n}$，因此样本方差就只剩下$n-1$个自由度了。

我们使用python画出不同自由度$n$下的密度函数图：

# 使用scipy计算pdf画图(非自定义函数)
from scipy.stats import chi2
x = np.linspace(0.01,30,10000)  
plt.plot(x, chi2.pdf(x,df=4),'r-', lw=5, alpha=0.6, label='chi2(4)',c='red')
plt.plot(x, chi2.pdf(x,df=6),'r-', lw=5, alpha=0.6, label='chi2(6)',c='blue')
plt.plot(x, chi2.pdf(x,df=10),'r-', lw=5, alpha=0.6, label='chi2(10)',c='orange')
plt.xlabel("X")
plt.ylabel("p (x)")
plt.legend()
plt.show()

# 使用卡方分布的定义演示卡方分布
from scipy.stats import norm
n = 10
chi2_list = []
for i in range(100000):
    x_i = norm.rvs(loc=0,scale=1,size=10)
    chi2_T = np.sum(np.square(x_i))
    chi2_list.append(chi2_T)
sns.distplot(chi2_list,color='blue')
plt.xlabel("x")
plt.ylabel("频率 p")
plt.title("n="+str(n))
plt.show()

该密度函数的图像是取非负值的偏态分布， 其数学期望等于自由度n， 方差等于 2 倍自由度即2n， 即 $E\left({\chi }^2\right)=n,\mathrm{Var}\left({\chi }^2\right)=2n$。

为什么${\chi }^2$分布很重要呢？我们在学习样本方差时，一直没有给出样本方差的抽样分布，是因为在不同的总体分布假设下，样本方差的抽样分布都是不一样的，而在正态分布总体假设下，样本方差经过变换可以与卡方分布产生关系，具体来说就是：

设 $x_1,x_2,\cdots ,x_n$ 是来自正态总体 $N\left(\mu ,{\sigma }^2\right)$ 的样本， 其样本均值和样本方差分别为
$$\bar{x}=\frac{1}{n}\sum _{i=1}^n{x}_i\text{ 和 }{s}^2=\frac{1}{n-1}\sum _{i=1}^n{\left(x_i-\bar{x}\right)}^2,$$
那么，$\frac{(n-1)s^2}{{\sigma }^2}\sim {\chi }^2(n-1)$。

我们使用python对这个结论进行验证：

# 假设总体为N（0，1），抽样的样本容量为n=11,抽样的次数为N
from scipy.stats import norm
def S2_Chi2(N):
    mu,sig = 0,1
    n = 11
    T_list = []
    for i in range(N):
        x_i = norm.rvs(loc=mu,scale=sig,size=n)  # 正态分布总体抽样
        T = (n-1)*np.var(x_i,ddof=1)/sig**2  # 构造卡方统计量
        T_list.append(T)
    sns.distplot(T_list,color='blue')
    plt.xlabel("x")
    plt.ylabel("频率 p")
    plt.title("chi2(10)")
    plt.show()

S2_Chi2(100000)

• （2.2）F统计量与F分布：

设随机变量 $X_1\sim {\chi }^2(m),X_2\sim {\chi }^2(n),X_1$ 与 $X_2$ 独立， 则称 $F=\frac{X_1/m}{X_2/n}$ 的分布是自由度为 $m$ 与 $n$ 的 $F$ 分布， 记为 $F\sim F(m,n)$， 其中 $m$ 称为分子自由度， $n$ 称为分母自由度。F分布的密度函数为：
$$\begin{array}{rl}{p}_F(y)& =\frac{\Gamma \left(\frac{m+n}{2}\right){\left(\frac{m}{n}\right)}^{\frac{m}{2}}{y}^{\frac{m}{2}-1}{\left(1+\frac{m}{n}y\right)}^{-\frac{m+n}{2}}}{\Gamma \left(\frac{m}{2}\right)\Gamma \left(\frac{n}{2}\right)}\cdot \end{array}$$

# 使用scipy与matplotlib绘制不同的m，n下的F分布的密度函数
from scipy.stats import f
x = np.linspace(0.01,5,10000)  
plt.plot(x, f.pdf(x,4,4000),'r-', lw=5, alpha=0.6, label='F(4,4000)',c='red')
plt.plot(x, f.pdf(x,4,10),'r-', lw=5, alpha=0.6, label='F(4,10)',c='blue')
plt.plot(x, f.pdf(x,4,4),'r-', lw=5, alpha=0.6, label='F(4,4)',c='orange')
plt.plot(x, f.pdf(x,4,1),'r-', lw=5, alpha=0.6, label='F(4,1)',c='yellow')
plt.xlabel("X")
plt.ylabel("p (x)")
plt.legend()
plt.show()

# 使用F统计量的定义演示：
from scipy.stats import norm
m,n = 4,4000
F_list = []
for i in range(100000):
    chi2_m_sample = np.sum(np.square(norm.rvs(loc=0,scale=1,size=m))) # 卡方m统计量
    chi2_n_sample = np.sum(np.square(norm.rvs(loc=0,scale=1,size=n))) # 卡方n统计量
    F_T = (chi2_m_sample/m) / (chi2_n_sample/n)  # # F（m，n）统计量
    F_list.append(F_T)
sns.distplot(F_list,color='blue')
plt.xlabel("x")
plt.ylabel("频率 p")
plt.title("F(4,4000)")
plt.show()

F分布的密度函数的图像是一个只取非负值的偏态分布。接下来，我们来看看之前学过的样本均值和样本方差与F分布的联系：

设 $x_1,x_2,\cdots ,x_m$ 是来自 $N\left({\mu }_1,{\sigma }_1^2\right)$ 的样本， $y_1,y_2,\cdots ,y_n$ 是来自 $N\left({\mu }_2,{\sigma }_2^2\right)$ 的样本， 且此两样本相互独立， 记：
$$s_x^2=\frac{1}{m-1}\sum _{i=1}^m{\left(x_i-\bar{x}\right)}^2,s_y^2=\frac{1}{n-1}\sum _{i=1}^n{\left(y_i-\bar{y}\right)}^2,$$
其中
$$\bar{x}=\frac{1}{m}\sum _{i=1}^m{x}_i,\bar{y}=\frac{1}{n}\sum _{i=1}^n{y}_i,$$
则有
$$F=\frac{s_x^2/{\sigma }_1^2}{s_y^2/{\sigma }_2^2}\sim F(m-1,n-1)$$
特别， 若 ${\sigma }_1^2={\sigma }_2^2$， 则 $F=s_x^2/s_y^2\sim F(m-1,n-1)$。

• （2.3）t分布及其统计量：

设随机变量 $X_1$ 与 $X_2$ 独立且 $X_1\sim N(0,1),X_2\sim {\chi }^2(n)$, 则称 $t=\frac{X_1}{\sqrt{X_2/n}}$ 的分布为自由度为 $n$ 的 $t$ 分布， 记为 $t\sim t(n)$。

下面，我们使用scipy绘制t分布的密度函数：

# 使用scipy与matplotlib绘制不同的n下的t分布的密度函数
from scipy.stats import t
from scipy.stats import norm 
x = np.linspace(-6,6,10000)  
plt.plot(x, t.pdf(x,4),'--', lw=5, alpha=0.6, label='t (4)',c='red')
plt.plot(x, norm.pdf(x,loc=0,scale=1),'r-', lw=5, alpha=0.6, label='N (0,1)',c='yellow')
plt.plot(x, t.pdf(x,100),'--', lw=5, alpha=0.6, label='t (100)',c='blue')
plt.xlabel("X")
plt.ylabel("p (x)")
plt.legend()
plt.show()

可以看到：当自由度较大 $($ 如 $n⩾30)$ 时, $t$ 分布可以用 $N(0,1)$ 分布近似（图中$t(100)$几乎与$N(0,1)$重合）。下面，我们使用样本均值和样本方差构造t统计量：

设 $x_1,x_2,\cdots ,x_n$ 是来自正态分布 $N\left(\mu ,{\sigma }^2\right)$ 的一个样本， $\bar{x}$ 与 $s^2$ 分别是该样本的样本均值与样本方差， 则有
$$t=\frac{\sqrt{n}(\bar{x}-\mu )}{s}\sim t(n-1)$$

这个统计量的由来十分具有故事性：在1980年以前，统计学主要的工作时解决社会统计，如人口问题等，后来生物统计问题异军突起，这类问题的主要特点是：数据量一般较大，所用的方法大多数是以中心极限定理为依据的，因此由于中心极限定理总是会联系到正态分布，正态分布在那时候大行其道。皮尔逊认为：正态分布是上帝赐予人类唯一正确的分布类型。到了20世纪初期，越来越多农业、工业试验的统计数据，这些统计数据的特点是：数据量一般不大，没办法使用中心极限定理近似。1898年，酿酒化学技师戈塞特平时的工作接触的数据量很小，只有几个，他通过大量的实验数据发现：$t=\frac{\sqrt{n}(\bar{x}-\mu )}{s}\sim t(n-1)$ 的分布与标准正态分布$N(0,1)$并不一致，但是由于数学和统计学功底不够，并不能解决这个问题。于是，他到皮尔逊那里学习，着重研究少量数据的统计分析，在1908年终于以Student为笔名发表了有关论文，并提出：设 $x_1,x_2,\cdots ,x_n$ 是来自正态分布 $N\left(\mu ,{\sigma }^2\right)$ 的一个样本， $\bar{x}$ 与 $s^2$ 分别是该样本的样本均值与样本方差， 则有：$t=\frac{\sqrt{n}(\bar{x}-\mu )}{s}\sim t(n-1)$。

$t$ 分布的发现在统计学史上具有划时代的意义, 因为$t$分布的出现打破了正态分布一统天下的局面，$t$分布 开创了小样本统计推断的新纪元， 小样本统计分析由此引起了广大统计科研工作者的重视。事实上，戈塞特的证明存在着漏洞， 费希尔（Fisher）注意到这个问题并于 1922 年给出了此问题的完整证明。

from scipy.stats import norm 
from scipy.stats import t 
t_list = []
for i in range(300000):
    mu,sigma2 = 0,1
    x_i = norm.rvs(loc=mu, scale=sigma2, size=5)
    x_mean = np.mean(x_i)
    x_s = np.std(x_i,ddof=1)
    t_T = np.sqrt(4)*(x_mean-mu) / x_s
    t_list.append(t_T)
sns.distplot(t_list,color='blue',label='t (4)')
x = np.linspace(-6,6,10000)  
plt.plot(x, norm.pdf(x,loc=0,scale=1),'r-', lw=5, alpha=0.6, label='N (0,1)',c='yellow')
plt.xlabel("x")
plt.ylabel("频率 p")
plt.title("t (4)")
plt.xlim(-6,6)
plt.legend()
plt.show()

四、参数估计之点估计的概念

在之前的学习中，我们主要学习了统计量以及统计量的抽样分布，回到最本质的问题，为什么需要学习统计量呢？引入统计量的目的就是想使用样本信息取推断总体的信息，而统计量只含有样本的信息，不含未知的总体参数。如何对未知的总体进行推断呢？在之前的讨论中，我们知道总体是一个随机变量，而描述随机变量最完整的莫过于随机变量的分布，而在实际中我们感兴趣的多是分布中未知参数有关，如：全中国人的身高服从均值为$\mu$、方差为${\sigma }^2$的正态分布$N(\mu ,{\sigma }^2)$，而总体中的均值$\mu$、方差${\sigma }^2$都是未知的，我们需要在总体中抽样，去估计总体中的未知参数均值$\mu$和方差${\sigma }^2$，因此总结下：

估计：通过样本统计量对总体分布的未知参数进行估计。

估计的方法有点估计与区间估计，点估计希望使用一个数估计总体中的位置参数，如$\mu =0$就是指使用一个数0去估计总体中的参数$\mu$，而0是由抽样后计算某个样本统计量得来的。换句话说，区间估计指的是使用一个区间估计总体中的参数，区间估计解决了点估计无法评价估计的精度的问题，这点我们在后面详细看，我们先学习点估计。

设 $x_1,x_2,\cdots ,x_n$ 是来自总体的一个样本， 用于估计未知参数 $\theta$ 的统计量 $\hat{\theta }=\hat{\theta }\left(x_1,x_2,\cdots ,x_n\right)$ 称为 $\theta$ 的估计量， 或称为 $\theta$ 的点估计， 简称估计。

点估计的方法有很多，接下来主要学习：矩估计和极大似然估计。事实上，还有很多其他的方法进行点估计，如：贝叶斯估计等。

五、参数估计之点估计的方法：矩估计

（1）总体矩和中心矩
在学习矩估计的方法前，我们先来看看什么是矩的概念。对于一个随机变量X来说，随机变量X的矩可以分为原点矩和中心矩，具体来说：

设 $X$ 为随机变量， $k$ 为正整数。 如果以下的数学期望都存在， 则称
$${\mu }_k=E\left(X^k\right)$$
为 $X$ 的 $k$ 阶原点矩。 称
$${\nu }_k=E(X-E(X){)}^k$$
为 $X$ 的 $k$ 阶中心矩。

显然，数学期望是随机变量的1阶原点矩，方差是随机变量的2阶中心矩。随机变量的矩是随机变量的一类数字特征，随机变量的原点矩刻画了随机变量$X$偏离原点$(0,0)$的程度，而中心矩描述了随机变量$X$偏离“中心”的程度，可以使用数学期望和方差做类比。

与此同时，一类常见的统计量就是样本矩，具体来说：
设 $x_1,x_2,\cdots ,x_n$ 是样本， $k$ 为正整数， 则统计量
$$a_k=\frac{1}{n}\sum _{i=1}^n{x}_i^k$$
称为样本 $k$ 阶原点矩。 特别地， 样本一阶原点矩就是样本均值。 统计量
$$b_k=\frac{1}{n}\sum _{i=1}^n{\left(x_i-\bar{x}\right)}^k$$
称为样本 $k$ 阶中心矩。 特别地， 样本二阶中心矩就是样本方差。

下面，我们来看看在python如何计算样本矩：

# 假设总体是标准正态分布，求3阶原点矩和中心矩
from scipy.stats import norm
import numpy as np
x_i = norm.rvs(loc=0, scale=1, size=10000)
a3 = np.mean(np.power(x_i,3))
b3 = np.mean(np.power((x_i-np.mean(x_i)), 3))
print("3阶原点矩：",a3)
print("3阶中心矩：",b3)

3阶原点矩： -0.03721661726047597
3阶中心矩： -0.018036837246877747

（2）矩估计
1900年，皮尔逊提出了替换原理，具体来说就是：使用样本矩（样本原点矩和样本中心矩）替换总体矩（原点矩和中心矩），如：使用样本均值$\bar{x}$替换总体均值$E(X)$、使用样本方差$s^2$替换总体方差$Var(X)$。

当然，我们也可以拓展下替换原理就是矩估计，如：使用样本均值$\bar{x}$估计总体均值$E(X)$、使用样本方差$s^2$估计总体方差$Var(X)$、用事件的频率估计事件的概率、使用样本分位数估计总体分位数等。

为什么矩估计能够让人接受呢？原因就是格利纹科定理：使用经验分布函数替换总体分布。

【例子】假设总体服从指数分布，其密度函数为
$$p(x;\lambda )=\lambda {\mathrm{e}}^{-\lambda x},x⩾0$$
从该总体中抽样1000个样本，估计总体分布的参数$\lambda$。

解：

$k=1$， 由于 $E(X)=1/\lambda$，即 $\lambda =1/E(X)$，因此，$\lambda$的矩估计为：$\hat{\lambda }=\frac{1}{\bar{x}}$

# 假设真实值lambda = 5
from scipy.stats import expon
real_lmd = 5
x_i = np.random.exponential(scale=1/real_lmd, size=1000)
print("矩估计为：",1/np.mean(x_i))

矩估计为： 4.917129994515655

六、参数估计之点估计的方法：极大似然估计

极大似然估计的思想非常有意思，充分利用了样本的二重性，即：可以把样本想象成黑盒子，打开前样本是一个随机变量，打开后就是确定的值。具体来说：极大似然估计就是利用已知的样本结果信息，反推最有可能（最大概率）导致这些样本结果出现的模型参数值。下面举例子说明：

假设我们有一个罐子（总体），罐子里有黑白两色球，假设黑球为记为1，白球记为0，球的比例未知，假设黑球的比例为$p$，白球的比例时$1-p$。为了估计总体罐子的参数$p$，我们对总体进行抽样，采样的结果如下：
$$1,1,0,1,1,1,0,0,1,1$$
为了使用我们的样本信息估计总体参数$p$，我们计算这组样本出现的概率（似然函数）
$$\begin{array}{cc}P(p)& =p\times p\times (1-p)\times p\times p\times p\times (1-p)\times (1-p)\times p\times p\\ & =p^7\times (1-p{)}^3\end{array}$$
由于不同的$p$会导致样本发生的概率发生改变，可能是“冥冥之中自有天意”，有一股无形的力量迫使我们从总体中采样刚好采到这个样本，因此样本发生的概率应该是最大才对。接下来，我们需要把样本发生的概率最大化：
$$ma{x}_{p}P(p)=p^7\times (1-p{)}^3$$
由于直接求解$n$个连乘的式子是十分困难的，但是我们比较习惯求解加法的式子，因此可以使用对数把以上的式子$P(p)$简化为连加式子$ln(P(p))$（对数似然函数），即：
$$ln(P(p))=7ln(p)\times 3ln(1-p)$$
最大值的求解需要使用导数的知识，即导函数为0：
$$\begin{array}{cc}& \frac{dln(P)}{dp}=0\\ & \frac{7}{p}-\frac{3}{1-p}=0\\ & p=0.7\end{array}$$
求解的最大值对应的$p=0.7$就是我们所要估计的p。下面我们使用python求解：

# 使用sympy演示极大似然估计的案例
from sympy import *
p = Symbol('p')  #定义总体参数
P_p = p**7*(1-p)**3  # 定义似然函数
lnP_p = ln(P_p) # 化简为对数似然
d_ln_P = diff(lnP_p, p) # 求导函数
p_hat = solve(d_ln_P, p) # 导函数为0
print("p的极大似然估计为：",p_hat)

p的极大似然估计为： [7/10]

下面，我们来实践一个更加有意义的案例：对正态总体 $N\left(\mu ,{\sigma }^2\right),\theta =\left(\mu ,{\sigma }^2\right)$ 是二维参数， 设有样本 $x_1,x_2,\cdots ,x_n$，求总体参数$\mu$和${\sigma }^2$。
$$\begin{array}{rl}L\left(\mu ,{\sigma }^2\right)=& \prod _{i=1}^n\left(\frac{1}{\sqrt{2\pi }\sigma }\exp \left\{-\frac{{\left(x_i-\mu \right)}^2}{2{\sigma }^2}\right\}\right)={\left(2\pi {\sigma }^2\right)}^{-n/2}\exp \left\{-\frac{1}{2{\sigma }^2}\sum _{i=1}^n{\left(x_i-\mu \right)}^2\right\}\\ & \ln L\left(\mu ,{\sigma }^2\right)=-\frac{1}{2{\sigma }^2}\sum _{i=1}^n{\left(x_i-\mu \right)}^2-\frac{n}{2}\ln {\sigma }^2-\frac{n}{2}\ln (2\pi )\end{array}$$
为了求解对数似然函数的最大化，需要$lnL(\mu ,{\sigma }^2)$对$\mu$和${\sigma }^2$求偏导函数，然后求偏导函数的零点。
$$\begin{array}{c}\frac{\partial \ln L\left(\mu ,{\sigma }^2\right)}{\partial \mu }=\frac{1}{{\sigma }^2}\sum _{i=1}^n\left(x_i-\mu \right)=0\\ \frac{\partial \ln L\left(\mu ,{\sigma }^2\right)}{\partial {\sigma }^2}=\frac{1}{2{\sigma }^4}\sum _{i=1}^n{\left(x_i-\mu \right)}^2-\frac{n}{2{\sigma }^2}=0.\end{array}$$
解这个方程组，就可以得到：
$$\hat{\mu }=\frac{1}{n}\sum _{i=1}^n{x}_i=\bar{x}$$
和
$${\hat{\sigma }}^2=\frac{1}{n}\sum _{i=1}^n{\left(x_i-\bar{x}\right)}^2=s_n^2$$
我们可以惊讶地发现，正态分布的$\mu$和${\sigma }^2$的极大似然估计竟然是样本均值$\bar{x}$和有偏样本方差$s_n^2$，这就将我们的估计与统计量连接起来了。一般来说，估计的结果都与从该分布抽样的样本组成的样本统计量有关，如样本均值$\bar{x}$样本方差$s^2$ 等等。

【例子】指数分布通常用来描述事件之间的时间，假如在银行排队，每个客户排队的等待时间服从参数为$\lambda$的指数分布。指数分布的密度函数为：
$$p(x)=\{\begin{array}{r}\lambda e^{-\lambda x},x⩾0\\ 0,x<0\end{array}$$
参数$\lambda$代表发生率，如：$\lambda =1$代表平均每分钟有1个人完成服务，$\lambda =2$代表平均1分钟有两位客户完成服务，$\lambda =0.5$代表每2分钟有1位客户完成服务（也可以说是每分钟有0.5位客户完成服务。）

为了估计银行的服务速度，我们随机在银行的客户进行抽样，询问本次服务的时间长度，得到5位客户的回答：$x_1=2,x_2=3,x_3=0.5,x_4=5,x_5=2$，下面我们使用极大似然估计去估计总体分布（指数分布）中的参数$\lambda$。

解：

似然函数：
$$L(\lambda |x_1,x_2,...x_5)={\lambda }^{\mathrm{n}}\left[e^{-\lambda \left(x_1+x_2+\dots +x_5\right)}\right]$$
对数似然函数为：
$$ln(L(\lambda |x_1,x_2,...,x_5))=5\log (\lambda )-\lambda \left(x_1+x_2+\dots +x_5\right)$$
对数似然函求导：
$$\frac{dln(L(\lambda |x_1,x_2,...,x_5))}{d\lambda }=n\frac{1}{\lambda }-\left(x_1+x_2+\dots +x_5\right)$$
导函数为0：
$$\begin{array}{cc}& n\frac{1}{\lambda }-\left(x_1+x_2+\dots +x_5\right)=0\\ & \lambda =\frac{n}{\left(x_1+x_2+\dots +x_5\right)}=\frac{1}{\bar{x}}\end{array}$$

from sympy.abc import lamda
x_1,x_2,x_3,x_4,x_5 = symbols('x_1:6') # 定义多个样本变量
x_1,x_2,x_3,x_4,x_5 = 2, 3, 0.5, 5, 2
f_lmd = lamda*E**(-lamda*x_1) * lamda*E**(-lamda*x_2) * lamda*E**(-lamda*x_3) * lamda*E**(-lamda*x_4) * lamda*E**(-lamda*x_5) # 定义似然函数
ln_f_lmd = ln(f_lmd) # 定义对数似然函数
d_ln_f = diff(ln_f_lmd, lamda) # 求导
lmd_hat = solve(d_ln_f, lamda) # 导数为0
print("指数分布参数lamda的极大似然估计值为：",lmd_hat)
print("指数分布参数lamda的极大似然公式求解为：", 5 / (x_1+x_2+x_3+x_4+x_5))

指数分布参数lamda的极大似然估计值为： [0.400000000000000, zoo]
指数分布参数lamda的极大似然公式求解为： 0.4

因此，$\lambda =0.4$表示银行队伍每分钟有0.4位客户完成服务，也就是说银行队伍2.5分钟有1人完成服务。