GitModel-动手学数理统计_04（python）

1 动手学数理统计_04

github 上pdf版本及ipynb版本：https://github.com/cx-333/Math-Modeling

1.11 假设检验之似然比检验与Bootstrap方法

• 似然比检验：

  定义（广义似然比）：设 $x_1,\cdots ,x_n$ 为来自密度函数为 $p(x;\theta ),\theta \in \Theta$ 的样本，考虑检验问题
$$H_0:\theta \in {\Theta }_{0}vs{H}_1:\theta \in {\Theta }_1=\Theta -{\Theta }_0$$
令
$$\Lambda \left(x_1,\cdots ,x_n\right)=\frac{{\sup }_{\theta \in \Theta }p\left(x_1,\cdots ,x_n;\theta \right)}{{\sup }_{\theta \in {\Theta }_0}p\left(x_1,\cdots ,x_n;\theta \right)}\sup 表示最小上界$$
那么称它为假设检验问题的广义似然比。

  由上式可以看出，广义似然比定义下的分子分母都有一个上确界（$\sup$）的符号，仔细观察，分子就相当于在全参数空间下取联合概率密度的最大值，分母相当于在原假设参数空间下取联合概率密度的 最大值，所以这个比值就是两个极大似然估计的比值。直观来看，如果原假设是正确的，那么参数 应该会落在原假设的参数空间内，换句话说，分子的最大值对应的参数应该落在 ${\Theta }_0$ ，所以这个 比值就不会太大。但是，反过来说，如果原假设应该被拒绝，那么参数就有很大可能落在拒绝域， 那么全参数空间的最大值就会在 $\theta \in {\Theta }_1$ 中取到，那么这个时候比值就会变大。所以可以看出 来，拒绝域顺理成章的应该设置为
W = { Λ ( x 1 , ⋯   , x n ) ≥ c } W=\left\{\Lambda\left(x_{1}, \cdots, x_{n}\right) \geq c\right\} W={Λ(x1​,⋯,xn​)≥c}
其中临界值 $c$ 要满足 $P_{\theta }\left(\Lambda \left(x_1,\cdots ,x_n\right)\ge c\right)\le \alpha ,\forall \theta \in {\Theta }_0$.

💡 $2ln\Lambda (\left(x_1,\cdots ,x_n\right))$服从${\chi }^2$分布，自由度为独立参数的个数（需要检验参数的维度）

🔥例子:假设观察某种疾病的发生情况: $n=100$ 人中发生了 $k=10$ 个事件。假定数据服从二项分布, 理论已知人群中每个人发生该事件的概率为 ${\pi }_0=0.2$ 。试对该假设做似然比假设检验（显著性水平$\alpha =0.05$）?

🦊解：

1. 作出假设
   $$H_0:\pi ={\pi }_0=0.2;H_1:\pi \ne {\pi }_0=0.2$$
2. 写出参数空间
   Θ 0 = { π 0 } , Θ = { π , π ∈ R } \Theta_{0}=\{\pi_{0}\}, \quad \Theta = \{\pi, \pi \in \textbf{R}\} Θ0​={π0​},Θ={π,π∈R}
3. 计算全参数空间（$\Theta$）下的极大似然估计
   $$\begin{array}{rl}& f(x)=C_n^k{p}^k(1-p{)}^{n-k}\\ & lnf(x)=ln{C}_n^k+ln{p}^k+ln(1-p{)}^{n-k}\\ & 令\frac{d}{dp}lnf(x)=0\\ & 解得\\ & p=\frac{k}{n}\end{array}$$
4. 在$H_0$假设下，$p=0.2$，计算$2ln\frac{{\sup }_{\theta \in \Theta }p(x_1,x_2,\cdots ,x_n)}{{\sup }_{\theta \in {\Theta }_0}p(x_1,x_2,\cdots ,x_n)}$，为
   $$check=2ln\frac{k/n}{0.2}$$
5. 因为
   $$2ln\frac{{\sup }_{\theta \in \Theta }p(x_1,x_2,\cdots ,x_n)}{{\sup }_{\theta \in {\Theta }_0}p(x_1,x_2,\cdots ,x_n)}\sim {\chi }^2(1)$$
6. 因此，拒绝域为
   $$check\ge {\chi }_{\alpha /2}(1)或check\le {\chi }_{1-\alpha /2}(1)$$
7. 判断$check$值是否在拒绝域，若在，则拒绝$H_0$，否则接受$H_0$.

python代码（求解上题）

import numpy as np 
from scipy.stats import chi2 

p0 = 0.2
n = 100
k = 10
alpha = 0.05
check = 2 * np.log((k/n)/p0)
left = chi2.ppf(df=1, q=alpha/2)
right = chi2.ppf(df=1, q=(1-alpha/2))
print("似然比统计量为:{}".format(check))
print("拒绝域为(0, {})或({}, oo)".format(left, right))
if check <= left or right >= right:
    print("拒绝假设H0, 即每个人发生该事件得概率不为0.2")

似然比统计量为:-1.3862943611198906
拒绝域为(0, 0.0009820691171752555)或(5.023886187314888, oo)
拒绝假设H0, 即每个人发生该事件得概率不为0.2

• Bootstrap方法：在上述方法中，如果分布的自由度难以确定，这个方法将难以进行下去。换句话说，当碰到某个统计量的分布难以确定或者未知的时候如何做假设检验呢？Bootstrap方法就是在这一背景下产生的。设总体得分布$F$未知，但已经有一个容量为$n$得来自分布$F$的数据样本，自这一样本按放回抽样的方法抽取一个容量为$n$的样本，这种样本称为bootstrap样本或自助样本，相继地、独立地自原始样本中抽取很多个bootstrap样本，利用这些样本对总体$F$进行统计推断，这种方法称为非参数bootstrap方法。

bootstrap置信区间：设$X=X_1,X_2,\cdots ,X_n$是来自总体$F$容量为$n$的样本，$x=x_1,x_2,\cdots ,x_n$是一个已知的样本值，$F$中含有未知参数$\theta$，$\hat{\theta }=\hat{\theta }(X_1,X_2,\cdots ,X_n)$是$\theta$的估计量，$\theta$的置信水平为$1-\alpha$（显著性水平为$\alpha$）的置信区间为:相继地、独立地从样本$x=x_1,x_2,\cdots ,x_n$中抽出B个容量为n的bootstrap样本，对于每个样本求出的$\theta$的bootstrap估计：$\widehat{{\theta }_1},\widehat{{\theta }_2},\cdots ,\widehat{{\theta }_n}$,将他们从小到大排序：$$\widehat{{\theta }_{(1)}}<\widehat{{\theta }_{(2)}}<\cdots <\widehat{{\theta }_{(n)}}$$ 置信区间取$$(\hat{\theta }(k_1),\hat{\theta }(k_2))$$其中，$k_1=[B\frac{\alpha }{2}],k_2=[B(1-\frac{\alpha }{2})],([\cdot ]表示取整)$

🔥例子：某工厂生产以发光产品，发光产品的发光时长服从正态分布 $N\left(\mu ,{\sigma }^2\right)$ ，产品的发光时长设定均值为 $250\mathrm{h}$ 。现在从一批产品中抽取 $10$ 个产品，测得发光时长为（单位 为: $h$) :

$$248.8,249.2,250.7,251.2,248.0,253.0,248.9,250.2,251.2,249.2$$

问该厂的发光产品是否符合要求（显著性水平$\alpha =0.05$）?

🦊解：该问题为左边检验，单侧置信区间的形式应为$(a,+oo)$

1. 作出假设
   $$H_0:\mu \ge {\mu }_0=250;H_1:\mu <{\mu }_0=250$$
2. 确定检验统计量$\overline{x}$
3. 进行$B$次bootstrap采样并计算参数${\mu }_i,(i=1,2,\cdots ,B)$后从小到大排序${\mu }_{(1)}<{\mu }_{(2)}<\cdots <{\mu }_{(B)}$
4. 取$\alpha /2,1-\alpha /2$分位点，得到置信区间$(\hat{\theta }(k_1),\hat{\theta }(k_2))（这是双边检验的置信区间，这里须修改为单边）$, $k_1=[B\frac{\alpha }{2}],k_2=[B(1-\frac{\alpha }{2})],([\cdot ]表示取整)$
5. 判断检验统计量是否在置信区间内，若在接受$h_0$，否则，拒绝$H_0$.

python代码(计算上题)

import numpy as np

# 例题为左边检验，单侧置信区间为（k1, oo）,拒绝域为(-oo, k1]

mu0 = 250
alpha = 0.05
B = 1000
x = [248.8, 249.2, 250.7, 251.2, 248.0, 253.0, 248.9, 250.2, 251.2, 249.2]
x_mean= np.mean(x)
params = []
for i in range(B):
    x_resample = np.random.choice(x, len(x), replace=True)
    params.append(np.mean(x_resample))

params = np.sort(params)

k1 = int(B*alpha)
left = params[k1-1]
# right = np.percentile(params, (1-alpha)*100-1)
print("检验统计量为:{}".format(x_mean))
print("拒绝域为:(-oo, {}]".format(left))
if x_mean > left:
    print("接受H0, 该厂的发光产品符合要求")
else:
    print("拒绝H0, 该厂的发光产品不符合要求")

检验统计量为:250.03999999999996
拒绝域为:(-oo, 249.32000000000002]
接受H0, 该厂的发光产品符合要求

$k_1=[B\frac{\alpha }{2}],k_2=[B(1-\frac{\alpha }{2})],([\cdot ]表示取整)$

🔥例子：某工厂生产以发光产品，发光产品的发光时长服从正态分布 $N\left(\mu ,{\sigma }^2\right)$ ，产品的发光时长设定均值为 $250\mathrm{h}$ 。现在从一批产品中抽取 $10$ 个产品，测得发光时长为（单位 为: $h$) :

$$248.8,249.2,250.7,251.2,248.0,253.0,248.9,250.2,251.2,249.2$$

问该厂的发光产品是否符合要求（显著性水平$\alpha =0.05$）?

🦊解：该问题为左边检验，单侧置信区间的形式应为$(a,+oo)$

1. 作出假设
   $$H_0:\mu \ge {\mu }_0=250;H_1:\mu <{\mu }_0=250$$
2. 确定检验统计量$\overline{x}$
3. 进行$B$次bootstrap采样并计算参数${\mu }_i,(i=1,2,\cdots ,B)$后从小到大排序${\mu }_{(1)}<{\mu }_{(2)}<\cdots <{\mu }_{(B)}$
4. 取$\alpha /2,1-\alpha /2$分位点，得到置信区间$(\hat{\theta }(k_1),\hat{\theta }(k_2))（这是双边检验的置信区间，这里须修改为单边）$, $k_1=[B\frac{\alpha }{2}],k_2=[B(1-\frac{\alpha }{2})],([\cdot ]表示取整)$
5. 判断检验统计量是否在置信区间内，若在接受$h_0$，否则，拒绝$H_0$.

python代码(计算上题)

import numpy as np

# 例题为左边检验，单侧置信区间为（k1, oo）,拒绝域为(-oo, k1]

mu0 = 250
alpha = 0.05
B = 1000
x = [248.8, 249.2, 250.7, 251.2, 248.0, 253.0, 248.9, 250.2, 251.2, 249.2]
x_mean= np.mean(x)
params = []
for i in range(B):
    x_resample = np.random.choice(x, len(x), replace=True)
    params.append(np.mean(x_resample))

params = np.sort(params)

k1 = int(B*alpha)
left = params[k1-1]
# right = np.percentile(params, (1-alpha)*100-1)
print("检验统计量为:{}".format(x_mean))
print("拒绝域为:(-oo, {}]".format(left))
if x_mean > left:
    print("接受H0, 该厂的发光产品符合要求")
else:
    print("拒绝H0, 该厂的发光产品不符合要求")

检验统计量为:250.03999999999996
拒绝域为:(-oo, 249.32000000000002]
接受H0, 该厂的发光产品符合要求