区间估计

区间估计经常用于质量控制领域来检测生产过程是否正常运行或者在“控制之中” ，也可以用来监控互联网领域各类数据指标是否在正常区间。本篇我们将介绍总体均值的区间估计和总体比例的区间估计。

区间估计的Python函数实现如下：

import numpy as np
import scipy.stats
from scipy import stats as sts


def mean_interval(mean=None, sigma=None,std=None,n=None,confidence_coef=0.95):
    """
    mean:样本均值
    sigma: 总体标准差
    std: 样本标准差
    n:   样本量
    confidence_coefficient：置信系数
    confidence_level:置信水平 置信度
    alpha:显著性水平
    功能：构建总体均值的置信区间
    """
    alpha = 1 - confidence_coef
    z_score = scipy.stats.norm.isf(alpha / 2)            # z分布临界值
    t_score = scipy.stats.t.isf(alpha / 2, df = (n-1) )  # t分布临界值
   
    if n >= 30: 
        if sigma != None:
            me = z_score * sigma / np.sqrt(n)
            print("大样本，总体 sigma 已知：z_score:",z_score)
        elif sigma == None:
            me = z_score * std / np.sqrt(n)
            print("大样本，总体 sigma 未知 z_score",z_score)
        lower_limit = mean - me
        upper_limit = mean + me
    if n < 30 :
        if sigma != None:
            me = z_score * sigma / np.sqrt(n)
            print("小样本，总体 sigma 已知 z_score * sigma / np.sqrt(n) \n z_score = ",z_score)
        elif sigma == None:
            me = t_score * std / np.sqrt(n)
            print("小样本，总体 sigma 未知 t_score * std / np.sqrt(n) \n t_score = ",t_score)
            
        print("t_score:",t_score)
        lower_limit = mean - me
        upper_limit = mean + me
        print("总体样本区间:[{0:.1f},{1:.1f}]".format(lower_limit, upper_limit))
    
    return (round(lower_limit, 1), round(upper_limit, 1))

• 网站独立访问量（UV）的区间估计示例

某网站流量UV数据如下[52,44,55,44,45,59,50,54,62,46,54,42,60,62,43,42,48,55,57,56]，我们研究一下该网站的总体流量uv均值：

data = np.array([52, 44, 55, 44, 45, 59, 50, 54, 62, 46, 54, 42, 60, 62, 43, 42, 48, 55, 57, 56])

# 计算均值
x_bar = data.mean()
print("均值：{0:.1f}".format(x_bar))
# 51.5

# 样本标准差
x_std = sts.tstd(data, ddof=1)  # ddof=1时,分母为n-1;ddof=0时,分母为n
print("样本标准差：{0:.1f}".format(x_std))
# 6.840283158189472

mean_interval(mean=x_bar, sigma=None, std=x_std, n=20, confidence_coef=0.95)

运行上面程序，输出结果为：

均值：51.5
样本标准差：6.8
小样本，总体 sigma 未知 t_score * std / np.sqrt(n) 
 t_score =  2.093024054408263
t_score: 2.093024054408263
流量区间:[48.3,54.7]

于是我们有95%的把握，该网站的流量uv介于 [48.3, 54.7]之间。接下去我们计算总体比例的区间估计，假定上述UV数据占据了总访问量的60%，即0.6，计算区间估计的函数如下：

# 总体比例的区间估计
def proportion_interval(p=None, n=None, confidence_coef=0.95):
    """
    p: 样本比例
    n: 样本量
    confidence_coef: 置信系数
    功能：构建总体比例的置信区间
    """
    alpha = 1 - confidence_coef
    z_score = scipy.stats.norm.isf(alpha / 2)  # z分布临界值

    me = z_score * np.sqrt((p * (1 - p)) / n)
    lower_limit = p - me
    upper_limit = p + me
    print("总体比例区间:[{0:.1f},{1:.1f}]".format(lower_limit, upper_limit))

    return (round(lower_limit, 3), round(upper_limit, 3))

调用程序与输出结果如下：

proportion_interval(p=0.6, n=20, confidence_coef=0.95)
# 总体比例区间:[0.4,0.8]

假设检验

统计领域要研究的是数据之间的区别和联系，也就是差异性和相关性分析。相关性分析在前面的描述性统计有叙述，这里重点关注数据的差异性。

比较两个数之间的大小，要么前后两者求差，要么求比。差值大于零说明前者大于后者，比值大于1说明分子大于分母。比较两组数据的差异性，原理和上面类似。从简单的开始，先比较一组数和一个给定数的差异，即单个总体的差异性分析。

• 单个总体的假设检验

常见单个总体差异性的假设检验分为均值、比例和方差：

A. 一个总体均值的假设检验（指定值和样本均值），

就是检验指定值与样本均值的差异，如下面食盐检验，按是否已知可以分2种情况：

1、均值已知情况：检验

假设我们要检验一批厂家生产的红糖是否够标重。监督部门称了50包标重500g的红糖，均值是498.35g，少于所标的500g。对于厂家生产的这批红糖平均起来是否够份量，需要统计检验。

分析过程：由于厂家声称每袋500g,因此原假设为总体均值等于500g(被怀疑对象总是放在零假设)。而且由于样本均值少于500g(这是怀疑的根据)，把备择假设设定为总体均值少于500g (上面这种备选假设为单向不等式的检验称为单侧检验，而备选假设为不等号“”的称为双侧检验，后面会解释)。于是原假设和备选假设如下：

H0：红糖均值等于500g

H1：红糖重量均值少于500g

import statsmodels.stats.weightstats

data = [493.01, 498.83, 494.16, 500.39, 497.63, 499.72, 493.41, 498.97, 501.94, 503.45, 497.47, 494.19, 500.99, 495.81,
        499.63, 494.91, 498.90, 502.43, 491.34, 497.50, 505.95, 496.56, 501.66, 492.02, 497.68, 493.48, 505.40, 499.21,
        505.84, 499.41, 505.65, 500.51, 489.53, 496.55, 492.26, 498.91, 496.65, 496.38, 497.16, 498.91, 490.98, 499.97,
        501.21, 502.85, 494.35, 502.96, 506.21, 497.66, 504.66, 492.11]

z, pval = statsmodels.stats.weightstats.ztest(data, value=500, alternative='smaller')

# 'two-sided': 样本均值与给定的总体均值不同
# 'larger' :   样本均值小于给定总体均值
# 'smaller' :  样本均值大于给定总体均值
print("z值={0:.4f},p值={1:.4f}".format(z, pval))
# z值=-2.6962,p值=0.0035

结论：选择显著性水平 0.05 的话，P=0.0035 < 0.05, 故应该拒绝原假设。具体来说就是该结果倾向于支持平均重量小于500g的备则假设。
2、均值未知情况：检验

假设我们检验汽车实际排放是否低于其声称的排放标准。汽车厂商声称其发动机排放标准的一个指标平均低于20个单位。在抽查了10台发动机之后,得到下面的排放数据:17.0 21.7 17.9 22.9 20.7 22.4 17.3 21.8 24.2 25.4该样本均值为21.13.究竟能否由此认为该指标均值超过20?

分析过程：由于厂家声称指标平均低于20个单位,因此原假设为总体均值等于20个单位(被怀疑对象总是放在零假设)。而且由于样本均值大于20(这是怀疑的根据)，把备择假设设定为总体均值大于20个单位

H0：汽车实际排放总体均值等于20个单位

H1：汽车实际排放总体均值大于20个单位

data = [17.0, 21.7, 17.9, 22.9, 20.7, 22.4, 17.3, 21.8, 24.2, 25.4]

t, pval = scipy.stats.ttest_1samp(a = data, popmean=20,alternative = 'greater')
# 说明
# a  为给定的样本数据
# popmean 为给定的总体均值
# alternative 定义备择假设。以下选项可用(默认为“two-sided”)：
# ‘two-sided’：样本均值与给定的总体均值(popmean)不同
# ‘less’：样本均值小于给定总体均值(popmean)
# ‘greater’：样本均值大于给定总体均值(popmean)

print("t值={0:.4f},p值={1:.4f}".format(t, pval))
# '''
# t值=1.2336,p值=0.1243
# '''

结论：选择显著性水平 0.01 的话，P=0.1243 > 0.05, 故无法拒绝原假设。具体来说就是该结果无法支持指标均值超过20的备则假设。

B. 一个总体比例的假设检验(指定比例和样本比例)

假设我们要检验高尔夫球场女性球员比例是否因促销活动而升高。某高尔夫球场去年打球 ?的人当中有20%是女性，为了增加女性球员的比例，该球场推出了一项促销活动来吸引更多的女性参加高尔夫运动，在活动实施了1个月后，球场的研究者想通过统计分析研究确定高尔夫球场的女性球员比例是否上升，收集到了400个随机样本，其中有100是女性。

分析过程：由于研究的是女性球员所占的比例是否上升，因此选择上侧检验比较合适,备择假设是比例大于20%

H0：女性球员所占比例小于等于20%

H1：女性球员所占比例大于20%

import statsmodels.stats.weightstats
import scipy.stats
import numpy as np
from statsmodels.stats.proportion import proportions_ztest
from scipy import stats as sts

count = 100
nobs = 400
p_0 = 0.2


p_bar = count/nobs
p_0 = 0.2
n = 400

# 执行单一样本比例检验 statsmodels.stats.proportion.proportions_ztest
#z_statistic, p_value = proportions_ztest(count, nobs, value=p_0, alternative='larger', prop_var=p_bar)
# 注：statsmodels.stats.proportion.proportions_ztest 的函数有几个问题：讲在第八节之后说明，感兴趣的读者请持续关注

def calc_z_score(p_bar, p_0, n):
    z = (p_bar - p_0) / (p_0 * (1 - p_0) / n)**0.5
    return z

z = calc_z_score(p_bar, p_0, n)
p = sts.norm.sf(z)

# 打印结果
print("z统计量：", z)
print("p值：", p)
#z统计量：2.4999999999999996
#p值： 0.006209665325776138

结论：选择显著性水平 0.05 的话，P=0.0062 < 0.05, 拒绝原假设。具体来说就是该结果支持特定的促销活动能够提升该球场女性运动员比例的备则假设。

C. 一个总体方差的假设检验（指定方差和样本方差）

假设我们检验公交车到站时间的方差是否比规定标准大。某市中心车站为规范化提升市民对于公交车到站时间的满意度，对于公交车的到站时间管理做了规定，标准是到站时间的方差不超过4。为了检验时间的到站时间的方差是否过大，随机抽取了24辆公交车的到站时间组成一个样本，得到的样本方差是，假设到站时间的总体分布符合正态分布，请分析总体方差是否过大。

分析过程：由于研究的是方差是否过大，因此选择上侧检验比较合适,备择假设是方差大于4，于是我们有了原假设和备用假设。

H0：公交到站的方差不超过4

H1：公交到站的方差大于4

def chi2test(sample_var, sample_num,sigma_square,side, alpha=0.05):
    '''
    参数：
    sample_var--样本方差
    sample_num--样本容量
    sigma_square--H0方差
    返回值：
    pval
    '''
    chi_square =((sample_num-1)*sample_var)/(sigma_square)
    p_value = None
    if side == 'two-sided':
        p = stats.chi2(df=sample_num-1).cdf(chi_square)
        p_value = 2*np.min([p, 1-p])
    elif side == 'less':
        p_value = stats.chi2(df=sample_num-1).cdf(chi_square)
    elif side == 'greater':
        p_value = stats.chi2(df=sample_num-1).sf(chi_square)
    return chi_square,p_value

chi_square,p_value = chi2test(sample_var = 4.9, sample_num = 24, sigma_square = 4,side='greater')

print("p值：", p_value)
# p值： 0.2092362676676498

结论：选择显著性水平 0.05 的话，P=0.2092 > 0.05, 无法拒绝原假设。具体来说就是该结果不支持方差变大的备则假设。

• 两个总体的假设检验

A. 两总体均值之差的假设检验(独立样本)

假设我们检验某药物在实验组的指标是否低于对照组。为检测某种药物对情绪的影响，对实验组的100名服药者和对照组的150名非服药者进行心理测试，得到相应的某指标。需要检验实验组指标的总体均值是否大于对照组的指标的总体均值。这里假定两个总体独立地服从正态分布。

分析过程：由于目标是检验实验组指标的总体均值是否大于对照组的指标的总体均值，因此选择上侧检验。于是我们有了原假设和备择假设

H0：实验组总计均值等于对照组均值

H1：实验组总体均值大于对照组均值

数据（drug.txt）如下，其中id=1为实验组，id=2为对照组

data = pd.read_table("./t-data/drug.txt",sep = ' ')
data.sample(5)

a = data[data['id']==1]['ah']
b = data[data['id']==2]['ah']

t, pval = scipy.stats.ttest_ind(a,b,alternative = 'greater')
print(t,pval)
# 0.9109168350628888 0.18161186154576608

结论：选择显著性水平 0.05 的话，p = 0.1816 > 0.05,无法拒绝H0，具体来说就是该结果无法支持实验组均值大于对照组的备则假设。

B. 两总体均值之差的假设检验(配对样本)

假设我们检验减肥前后的重量是否有显著性差异（是否有减肥效果）。这里有两列50对减肥数据（diet.txt）。其中一列数据(变量名before)是减肥前的重量，另一列(变量名after)是减肥后的重量(单位: 公斤)，人们希望比较50个人在减肥前和减肥后的重量。

分析过程：这里不能用前面的独立样本均值差的检验，这是因为两个样本并不独立。每一个人减肥后的重量都和自己减肥前的重量有关，但不同人之间却是独立的，所以应该用配对样本检验。同时，由于研究的是减肥前后的重量变化，期望减肥前的重量大于减肥后的重量，所以备择假设是期望减肥前的重量大于减肥后的重量

H0：期望减肥前的重量等于减肥后的重量

H1：期望减肥前的重量大于减肥后的重量

data = pd.read_table("./t-data/diet.txt",sep = ' ')
data.sample(5)

a = data['before']
b = data['after']
stats.ttest_rel(a, b,alternative = 'greater')
# Ttest_relResult(statistic=3.3550474801424173, pvalue=0.000769424325484219)

结论：选择显著性水平 0.05 的话，p = 0.0007 < 0.05,故应该拒绝原假设。具体来说就是该结果倾向支持减肥前后的重量之差大于零（即减肥前重量大于减肥后，也就是有减肥效果）的备则假设。

C. 两总体比例之差的假设检验

假设我们检验不同保险客户的索赔率是否存在差异。某保险公司抽取了单身与已婚客户的样本，记录了他们在一段数据内的索赔次数，计算了索赔率，现在需要检验两种保险客户的索赔率是否存在差异。

分析过程：由于目标比例是否有差异，因此选择比例之差的双侧检验

H0：两种客户索赔率不存在差异。

H1：两种客户索赔率存在差异。

def proportion_test(p1, p2, n1, n2, side='two-sided'):
    """
    参数：
    p1: 样本1的比例
    p2: 样本2的比例
    n1: 样本1的数量
    n2: 样本2的数量
    side: 假设检验的方向，可选'two-sided'（双侧检验，默认）, 'greater'（右侧检验）, 'less'（左侧检验）

    返回值：
    z_value: Z统计量的值
    p_value: 对应的p值
    """
    p = (p1 * n1 + p2 * n2) / (n1 + n2)
    se = np.sqrt(p * (1 - p) * (1 / n1 + 1 / n2))
    z_value = (p1 - p2) / se

    if side == 'two-sided':
        p_value = 2 * (1 - stats.norm.cdf(np.abs(z_value)))
    elif side == 'greater':
        p_value = 1 - stats.norm.cdf(z_value)
    elif side == 'less':
        p_value = stats.norm.cdf(z_value)
    else:
        raise ValueError("Invalid side value. Must be 'two-sided', 'greater', or 'less'.")

    return z_value, p_value


p1 = 0.14
p2 = 0.09
n1 = 250
n2 = 300

z_value, p_value = proportion_test(p1, p2, n1, n2, side='two-sided')
# 选择双侧检验 alternative = 'two-sided'

print("Z_value:", z_value)
print("p_value:", p_value)
# Z_value: 1.846189280616294
# p_value: 0.0648647268570739

结论：选择显著性水平 0.05 的话，p = 0.0648 > 0.05,故应该拒绝原假设。具体来说就是该结果倾向支持两种保险客户的索赔率存在差异的备则假设。

D. 两总体方差之比的假设检验

假设我们检验修完Python课程的学生是否比修完数据库课程的学生考CDA的成绩方差更大。某高校数据科学专业的学生，修完一门数据库课程的41名学生考CDA的方差，修完Python课程的31名学生考CDA的方差是，这些数据是否表明，修完数据库的学生要比修完Python的学生CDA成绩的方差更大？

分析过程：由于目标是希望修完Python的学生CDA成绩的方差更大，因此选择上侧检验。两总体方差之比用F检验,将方差较大的数据库课程的考试成绩视为总体1

H0：修完Python的学生CDN成绩方差更小

H1：修完Python的学生CDA成绩的方差更大

def f_test_by_s_square(n1, n2, s1_square,s2_square, side ='two-sided'):
    """
    参数
    n1 :样本1的数量
    n2 :样本2的数量
    s1_square:样本1的方差
    s2_square:样本2的方差
    # 
    # F_value :F统计量的值
    # p_value :对应的p值
    """
    F_value = s1_square/s2_square
    F = stats.f(dfn = n1-1, dfd = n2-1)
    if side=='two-sided':
        print("two-sided")
        p_value = 2*min(F.cdf(F_value), 1-F.cdf(F_value))
        return F_value,p_value
    elif  side=='greater':
        print("greater")
        p_value = 1-F.cdf(F_value)
        return F_value,p_value

f_statistic , p_value= f_test_by_s_square(n1=41, n2=31,s1_square=120,s2_square=80,side='greater')# 打印检验结果
# 选择上侧检验所以side='greater'
print("F statistic:", f_statistic)
print("p-value:", p_value)
# p-value：0.1256

结论：选择显著性水平 0.05 的话，p = 0.1256 > 0.05，故无法原假设。结果无法支持修完数据库的学生要比修完Python的学生CDA成绩的方差更大的备则假设。