推断统计学一般有两种方法，一是使用置信区间估算总体的参数，二是对总体参数的假设值进行决策。后者被称为假设检验，是我们这篇文章所要探讨的主题。
本文首先会介绍假设检验的逻辑，然后利用BRFSS数据探讨富人是否更瘦这个问题，除了使用t检验方法外，还采用了蒙特卡洛模拟的方法来实现假设检验。
假设检验的逻辑
假设检验的逻辑有些类似于反证法，这里有原假设和备择假设，它们之间是相互对立的，我们先在原假设成立的前提下进行计算，如果最终结果证明原假设不成立，我们就会接受备择假设，所以一般我们也会将想要证明的结论作为备择假设。
比如，在讨论富人是否比普通人更瘦的问题上，原假设是：富人和普通人一般胖瘦；而备择假设是：富人比普通人更瘦。然后我们通过实际的样本数据，证明在原假设成立的前提下：如果得到现有样本结果（甚至更极端的情况）的可能性微乎其微，以至于我们需要拒绝原假设而接受备择假设；如果可能性在我们能够接受的范围内，我们则无法拒绝原假设。
假设检验一般有下面几个步骤：
- 设置原假设和备择假设。
- 定义检验统计量。
- 在原假设的基础上，通过计算或模拟得到检验统计量和P值。
- 做判断：如果P值小于你能容忍的显著性水平时，则拒绝原假设。
这里的检验统计量是由样本数据计算而来的，如果检验统计量的值出现在拒绝域中，则拒绝原假设。而一般我们会计算检验统计量出现在极端情况下的概率，也就是P值，并与事先设置的显著性水平比较，如果P值更小，则意味着检验统计量处在拒绝域，需要拒绝原假设。
单个样本均值的假设检验
在开始假设检验前，我们需要先导入相应的Python模块和BRFSS数据，其中采用了scipy的stats统计模块，它提供了假设检验计算P值的函数。
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import brfss
from scipy import stats # 导入scipy的统计模块
%config InlineBackend.figure_format = 'retina'
df = brfss.ReadBrfss() # 读取数据
然后根据收入水平选取富人和普通人的bmi数据值。
df2 = df[['bmi', 'income']].dropna()
bmi = df2.bmi
bmi_rich = df2[df2.income == 8].bmi # 富人bmi数据
bmi_ord = df2[df2.income != 8].bmi # 普通人bmi数据
我们需要判断的问题是：富人会更瘦吗？为了演示单个样本的情况，这里只使用富人bmi数据和所有人总体bmi均值做比较。并假设总体的bmi均值等于样本bmi均值，由于样本量有40万，根据中心极限定理可以做这样的近似。
mu = np.mean(bmi) # 计算bmi均值
mu
28.18812531332513
在只使用富人bmi数据的单样本情况下，原假设为：富人bmi均值等于总体bmi均值28.188；备择假设为：富人bmi均值小于总体bmi均值28.188，这是一个单边检验。这里采用t检验，所以检验统计量为t统计量。
t检验一般用于变量满足或近似正态分布且总体标准差未知的情况下，特别是在小样本量的时候。但在大样本情况下，t检验和z检验是近似的，因此也可使用t检验。在scipy.stats模块中t检验的函数是`ttest_1samp()`，用于单个样本和总体均值比较的情况。其返回结果是t统计量和P值。
>```python
t_stats, p_value = stats.ttest_1samp(bmi_rich, mu)
print("p value is %.10f" % (p_value/2)) # 除以2是因为我们使用的是单边检验
p value is 0.0000000000
以上使用了所有富人的bmi数据，得到P值在目前的精度下为0，所以我们拒绝原假设，接受备择假设，即得出结论：富人更瘦。如果我们只使用部分bmi数据，则会得到完全不一样的结果。
t_stats, p_value = stats.ttest_1samp(bmi_rich[:500], mu) # 选择前500个数据进行检验
print("p value is %.10f" % (p_value/2))
p value is 0.2906513576
>```python
t_stats, p_value = stats.ttest_1samp(bmi_rich[:4000], mu) # 选择前4000个数据进行检验
print("p value is %.10f" % (p_value/2))
p value is 0.0000660545
当选择前500个数据进行t检验时，P值为0.29，而如果我们选择的显著性水平为0.1时，则无法拒绝原假设，但也不能证明原假设就成立。当选择前4000个数据进行t检验，P值为0.000066， 显著性水平为0.1时，则可以拒绝原假设。所以样本量越大，你得到的结果越明确。
蒙特卡洛模拟之 bootstrap 方法
以上直接采用了t检验，有时t检验的前提无法满足，我们还可以采用计算机模拟的方法，常用的是bootstrap方法。
bootstrap是一种自助的方法，即在原有的样本中进行可重复的重新抽样，得到更多的样本。因为我们先假设富人和普通人的bmi值没有差别，所以从富人的样本数据中进行重新抽样，其实是模拟了普通人的样本分布。然后将原富人样本bmi均值与多个新样本的bmi均值比较，从而计算出极端情况下的概率，即P值。
以下是实现bootstrap方法的代码。
def bootstrap_replicate_1d(data, func): # 进行一次重新抽样，并返回检验统计量
return func(np.random.choice(data, size=len(data)))
def draw_bs_reps(data, func, size=1):
bs_replicates = np.empty(size) # 初始一个空数组
for i in range(size): # 进行多次重新抽样
bs_replicates[i] = bootstrap_replicate_1d(data, func)
return bs_replicates # 返回多次抽样的检验统计量数组
def bootstrap_pvalue_1samp(data, pop_stats, func, size=1):
sample_stats = func(data) # 计算原有样本的检验统计量
translated_data = data - sample_stats + pop_stats # 数据平移
bs_replicates = draw_bs_reps(translated_data, func, size) # 重新抽样
p = np.sum( bs_replicates < sample_stats) / size # 计算抽样统计量小于原有统计量的概率
return p
同样，我们使用三组不同大小的数据得到了不同的结果，在大样本的情况下P值很小，需要拒绝原假设。
bootstrap_pvalue_1samp(bmi_rich, mu, np.mean, size=10000)
0.0
>```python
bootstrap_pvalue_1samp(bmi_rich[:500], mu, np.mean, size=10000)
0.71260000000000001
bootstrap_pvalue_1samp(bmi_rich[:4000], mu, np.mean, size=10000)
0.0001
### 两个独立样本均值的假设检验
现在我们讨论两个独立样本的情况，即使用富人bmi数据和普通人bmi数据。原假设为：富人和普通人的bmi均值相同；备择假设为：富人bmi均值比普通人小，仍是单边检验。同样，我们采用t检验，使用scipy.stats中的`ttest_ind()`函数，计算两个独立样本t检验的P值。
选择显著性水平为0.1时，当样本量只有500时，P值大于0.1，无法拒绝原假设。而当样本量为1000或全部时，P值小于0.1，拒绝原假设，证明了富人比普通人更瘦。
>```python
t_stats, p_value = stats.ttest_ind(bmi_rich, bmi_ord)
print("p value is %.10f" % (p_value/2))
p value is 0.0000000000
t_stats, p_value = stats.ttest_ind(bmi_rich[:500], bmi_ord[:500])
print("p value is %.10f" % (p_value/2))
p value is 0.1219871502
>```python
t_stats, p_value = stats.ttest_ind(bmi_rich[:1000], bmi_ord[:1000])
print("p value is %.10f" % (p_value/2))
p value is 0.0000846839
蒙特卡洛模拟之Permutation方法
同样，在两个独立样本的情况下，我们也可以采用计算机模拟的方法，这里介绍另一种Permutation的方法。
Permutation方法的逻辑是，在原假设成立的前提下，富人bmi和普通人bmi数据没有区别，所以可以将它们混合重新洗牌，然后重新分成两组数据，并计算这两组新数据均值的差值，重复多次这一操作。然后再将原有富人和普通人bmi均值之差与这些新组合的差值进行比较，计算极端情况下的概率，如果概率很小，则推翻了无差别的原假设。
以下是实现Permutation方法的代码。
def diff_of_means(data_1, data_2):
diff = np.mean(data_2) - np.mean(data_1) # 计算两组数据均值的差
return diff
def permutation_sample(data1, data2): # 产生新的分组数据
data = np.concatenate((data1, data2)) # 合并两组数据
permuted_data = np.random.permutation(data) # 对合并后的数据进行重新排列
perm_sample_1 = permuted_data[:len(data1)] # 分成新的两组数据
perm_sample_2 = permuted_data[len(data1):]
return perm_sample_1, perm_sample_2
def draw_perm_reps(data_1, data_2, func, size=1): # 进行多次重新分组的操作
perm_replicates = np.empty(size)
for i in range(size):
perm_sample_1, perm_sample_2 = permutation_sample(data_1, data_2)
perm_replicates[i] = func(perm_sample_1, perm_sample_2)
return perm_replicates
def permutation_pvalue(data_1, data_2, func, size=1): # 计算P值
empirical_test_stats = func(data_1, data_2)
perm_replicates = draw_perm_reps(data_1, data_2, func, size)
p = np.sum(perm_replicates > empirical_test_stats) / len(perm_replicates)
return p
同样，我们使用了大小不同的三组样本数据进行测试，得到与上述t分布一致的判断结果。
permutation_pvalue(bmi_rich, bmi_ord, diff_of_means, size=10000)
0.0
>```python
permutation_pvalue(bmi_rich[:500], bmi_ord[:500], diff_of_means, size=10000)
0.1142
permutation_pvalue(bmi_rich[:1000], bmi_ord[:1000], diff_of_means, size=10000)
0.0001
### 小结
至此，本系列的数据探索文章，使用美国行为风险因素监控BRFSS调研数据，围绕“富人会更瘦还是更胖” 这个问题，采用了多种统计学方法来寻找答案，它们分别是：
* 计算bmi均值的差
* 绘制统计图表
* 计算Cohen's d
* 计算bmi的置信区间
* 假设检验方法
由简到难，在解决问题的道路上，我们也将基础的统计方法在Python实现的情况下进行了整理和回顾，希望对大家有帮助。
---
[本文代码github](https://github.com/fishstar/Exploratory-Data-Analysis)
数据探索系列目录：
1. [开篇：数据分析流程](http://www.jianshu.com/p/dffdaf11bd4c)
2. [描述性统计分析](http://www.jianshu.com/p/8982ad63eb85)
3. [统计分布](http://www.jianshu.com/p/8a0479f55b21)
4. [参数估计](http://www.jianshu.com/p/7e556f17021a)
5. 假设检验（本文）
参考资料：
* [《Think Stats 2》](http://greenteapress.com/wp/think-stats-2e/)
* [《统计学》](https://book.douban.com/subject/4076214/)，William Mendenhall著
致谢：
最后还是非常感谢解密大数据社群的小伙伴们给我的支持和鼓励，让我们一起成长。
网友评论
原假设 和 备择假设 应该是完备且对立。文中所陈述明显有问题。