python 线性回归统计检验 p值_机器学习必备技能之“统计思维2.0”

weixin_39938724

902人浏览 · 2020-10-21 19:24:06

weixin_39938724 · 2020-10-21 19:24:06 发布

TUSHARE 金融与技术学习兴趣小组

编译整理 | 一只小绿怪兽

译者简介：北京第二外国语学院国际商务专业研二在读，目前在学习Python编程和量化投资相关知识。

在文章《机器学习必备技能之“统计思维1.0”》中，我们介绍了一些基本的统计工具，包括画图分析、定量分析、离散型随机变量和连续型随机变量四个部分，相信大家也已经掌握了这些技能，开始培养了用概率思考问题的统计思维。那么接下来，就可以着手深入研究数据集，解决实际问题了！

本文会拓展更多Python中的统计工具，并配合实例，向大家介绍在统计推断中两个非常关键的任务：

① 参数估计

② 假设检验

【工具】Python 3

【数据】tushare.pro 、Datacamp

【注】本文注重的是方法的讲解，请大家灵活掌握。

最优参数估计

我们发现，无论是离散型随机变量还是连续型随机变量，在公式中都有参数，选择的参数不同，结果也会不同。所以，最优参数估计的目的，就是找到最能够反映真实数据特征的参数。

以指数分布为例，它的概率密度函数中有一个参数λ，我们来看一下，如果给参数λ赋不同的值，它的分布函数对真实数据的拟合程度会有什么变化。

nohitter_times = np.array([843, 1613, 1101,  215,  684,  814,  278,  324,  161,  219,  545,
        715,  966,  624,   29,  450,  107,   20,   91, 1325,  124, 1468,
        104, 1309,  429,   62, 1878, 1104,  123,  251,   93,  188,  983,
        166,   96,  702,   23,  524,   26,  299,   59,   39,   12,    2,
        308, 1114,  813,  887,  645, 2088,   42, 2090,   11,  886, 1665,
       1084, 2900, 2432,  750, 4021, 1070, 1765, 1322,   26,  548, 1525,
         77, 2181, 2752,  127, 2147,  211,   41, 1575,  151,  479,  697,
        557, 2267,  542,  392,   73,  603,  233,  255,  528,  397, 1529,
       1023, 1194,  462,  583,   37,  943,  996,  480, 1497,  717,  224,
        219, 1531,  498,   44,  288,  267,  600,   52,  269, 1086,  386,
        176, 2199,  216,   54,  675, 1243,  463,  650,  171,  327,  110,
        774,  509,    8,  197,  136,   12, 1124,   64,  380,  811,  232,
        192,  731,  715,  226,  605,  539, 1491,  323,  240,  179,  702,
        156,   82, 1397,  354,  778,  603, 1001,  385,  986,  203,  149,
        576,  445,  180, 1403,  252,  675, 1351, 2983, 1568,   45,  899,
       3260, 1025,   31,  100, 2055, 4043,   79,  238, 3931, 2351,  595,
        110,  215,    0,  563,  206,  660,  242,  577,  179,  157,  192,
        192, 1848,  792, 1693,   55,  388,  225, 1134, 1172, 1555,   31,
       1582, 1044,  378, 1687, 2915,  280,  765, 2819,  511, 1521,  745,
       2491,  580, 2072, 6450,  578,  745, 1075, 1103, 1549, 1520,  138,
       1202,  296,  277,  351,  391,  950,  459,   62, 1056, 1128,  139,
        420,   87,   71,  814,  603, 1349,  162, 1027,  783,  326,  101,
        876,  381,  905,  156,  419,  239,  119,  129,  467])


def ecdf(data):
    """Compute ECDF for a one-dimensional array of measurements."""
    n = len(data)
    x = np.sort(data)
    y = np.arange(1, len(data)+1) / n

    return x, y

# 计算数据集nohitter_times的均值
mean = np.mean(nohitter_times)


# 模拟指数分布：参数取：均值
inter_nohitter_time = np.random.exponential(mean, 100000)
# 参数取：均值的一半
samples_half = np.random.exponential(mean/2, size=10000)
# 参数取：均值的两倍
samples_double = np.random.exponential(mean*2, size=10000)


# 用真实数据计算的 ECDF：x, y
x, y = ecdf(nohitter_times)

# 用估计值计算 ECDF: x_theor, y_theor
x_theor, y_theor = ecdf(inter_nohitter_time) 
x_half, y_half = ecdf(samples_half)
x_double, y_double = ecdf(samples_double)

# 画图：CDFs
plt.plot(x_theor, y_theor, label='mean')
plt.plot(x, y, marker='.', linestyle='none')

_ = plt.plot(x_half, y_half, label='mean_half')
_ = plt.plot(x_double, y_double, label='mean_double')

plt.margins(0.02)
plt.xlabel('Games between no-hitters')
plt.ylabel('CDF')

plt.legend()
plt.show()

很明显，当参数λ取样本均值的时候，对真实数据的拟合程度最好，可以作为最优参数使用。

寻找最优参数涉及到优化问题，模型不同，方法也不尽相同。本文会介绍一下大家比较熟悉的线性回归模型，以及它的最优参数估计方法：普通最小二乘法OLS。

关于线性回归的知识及应用，在文章《数据科学必备基础之线性回归》中已经进行了非常详细的介绍，这里就再简单回顾一下，重点是掌握最优化参数的思路。

一元线性回归模型可写为：y = α + βx，主要用于研究自变量x和因变量y之间的相关关系，其中β代表斜率，α代表截距，这里的α和β就是参数，它的图形是一条直线，如下所示。

我们知道，两点可以确定一条直线，但在真实的数据集中有很多点，这个时候该怎么确定直线的位置呢？

普通最小二乘法OLS(Ordinary Least Squares)提供了一个解决方案，思路就是找到一条直线，使得此直线距离所有这些点最近。这里的距离用残差平方和SSR(Sum of Squared Residuals)衡量。

残差(residual)是每个点到拟合直线的距离，如下图所示，残差平方和就是把所有的残差，先平方再加总。

这样，普通最小二乘法OLS就可以量化为，求最小化的残差平方和，由此得到的斜率β和截距α，就是最优参数估计值，数学表达式如下：

值得注意的是，在使用线性回归模型之前，对数据集进行探索是非常必要的，最直观的方法就是画图，这么做的目的是，判断数据集是否符合线性相关关系，否则就应该选择其他模型，或者做一些处理。

图片来自网络

比如，在上面四组图中，计算出的残差平方和是差不多的，但很显然除了set 1，线性回归模型并不是最佳选择。所以，在选择模型的时候，需要根据数据集的特征，慎重地做出决定，这就是提前做数据探索的意义，非常重要。

置信区间

统计推断依赖于概率，我们一直强调，用概率思考问题是非常重要的。

之前我们介绍过如何用某一个数据集的统计特征作为参数的估计值，然后根据这个估计值模拟出概率分布。但这里有一个问题，这个数据集能代表总体吗？如果换一组数据集，均值和标准差还会相同吗？答案大概率是不会。

所以，更好的一种选择是，我们随机模拟n次试验，然后用新生成的数据集计算统计值，这样计算出来的结果更具有代表性。

Bootstrap自举法是一个随机抽样的方法，它的原理是用重复抽样得到的数据，计算统计值。

第一步，我们可以直接调用np.random.choice()【1】函数对原数组进行随机抽样，效果如下。

import numpy as np

print(np.random.choice([1, 2, 3, 4, 5], size=5))
[3 3 3 3 2]

print(np.random.choice([1, 2, 3, 4, 5], size=3))
[1 3 1]

第二步，通过循环实现重复抽样。

import numpy as np
import tushare as ts
import matplotlib.pyplot as plt


ts.set_token('your token')
pro = ts.pro_api()

df = ts.pro_bar(ts_code='000001.SZ', adj='qfq', start_date='20190101', end_date='20190630')
df.sort_values('trade_date', inplace=True)
close_list = df['close'].tolist()

def ecdf(data):
    """Compute ECDF for a one-dimensional array of measurements."""
    n = len(data)
    x = np.sort(data)
    y = np.arange(1, len(data)+1) / n

    return x, y

for i in range(50):
    # 重复抽样
    bs_sample = np.random.choice(close_list, size=len(close_list))
    # 画图
    x, y = ecdf(bs_sample)
    _ = plt.plot(x, y, marker='.', linestyle='none', color='gray', alpha=0.1)

# 画图：原始数据
x, y = ecdf(close_list)
_ = plt.plot(x, y, marker='.')

plt.margins(0.02)
_ = plt.xlabel('000001.SZ_close_price')
_ = plt.ylabel('ECDF')
plt.show()

接下来，我们介绍一下置信区间(Confidence interval)的概念，同样涉及到概率思维。

在统计学中，一个概率样本的置信区间，是对这个样本的某个总体参数的区间估计。置信区间展现的是这个参数的真实值有一定概率落在测量结果的周围的程度，其给出的是被测量参数的测量值的可信程度，即前面所要求的“一个概率”，也叫置信水平。

也就是说，由于统计抽样存在误差，需要给计算出的统计值附上一个概率。

用同一个数据集close_list，我们可以试着计算一下收盘价均值的置信区间，即回答：有多大概率，收盘价的均值会在区间内。

首先，定义两个函数，实现bootstrap随机重复n次抽样，并计算出所需的统计值。

def bootstrap_replicate_1d(data, func):
    """
    随机抽样
    """
    return func(np.random.choice(data, size=len(data)))


def draw_bs_reps(data, func, size=1):
    """Draw bootstrap replicates."""

    bs_replicates = np.empty(size)

    for i in range(size):
        bs_replicates[i] = bootstrap_replicate_1d(data, func)

    return bs_replicates

随机重复抽样10000次，计算每次抽样的平均值，并画出概率分布图，可以发现，平均值是服从正态分布的。

bs_replicates = draw_bs_reps(close_list, np.mean, size=10000)
print(bs_replicates)

[12.17974576 12.11440678 12.15355932 ... 11.92813559 12.30347458
 12.2159322 ]

# 画图
_ = plt.hist(bs_replicates, bins=50, density=True)
_ = plt.xlabel('000001.SZ_close_price')
_ = plt.ylabel('PDF')
plt.show()

接下来，计算95%的置信水平下的置信区间，可以直接调用求分位数的函数np.percentile()。

# 计算95%的置信区间
conf_int = np.percentile(bs_replicates, [2.5, 97.5])
print('95% confidence interval =', str(conf_int))

95% confidence interval = [11.94101695 12.39077119]

前面我们介绍的是对数据集本身做随机重复抽样，且只涉及到一个平均值统计量，下面我们看一下如何用bootstrap自举法，来模拟包含两个参数的线性回归模型。

首先，定义一个函数，实现对线性回归模型的随机重复抽样。这里需要借助np.arange()【2】，生成一个有相同间隔的整数数组，为后面生成随机数做准备。

ts.set_token('your token')
pro = ts.pro_api()


df = ts.pro_bar(ts_code='000001.SZ', adj='qfq', start_date='20190101', end_date='20190630')
df.sort_values('trade_date', inplace=True)
# ----股票'000001.SZ'的收盘价数据
x = np.array(df['close'].tolist())


df_index = pro.index_daily(ts_code='399300.SZ', start_date='20190101', end_date='20190630')
df_index.sort_values('trade_date', inplace=True)
# ----沪深300指数的收盘价数据
y = np.array(df_index['close'].tolist())


def draw_bs_pairs_linreg(x, y, size=1):
    """Perform pairs bootstrap for linear regression."""

    inds = np.arange(0, len(x))

    bs_slope_reps = np.empty(size)
    bs_intercept_reps = np.empty(size)

    for i in range(size):
        bs_inds = np.random.choice(inds, size=len(inds))
        bs_x, bs_y = x[bs_inds], y[bs_inds]
        bs_slope_reps[i], bs_intercept_reps[i] = np.polyfit(bs_x, bs_y, 1)

    return bs_slope_reps, bs_intercept_reps

然后，调用函数计算在95%置信水平下，斜率的置信区间。

bs_slope_reps, bs_intercept_reps = draw_bs_pairs_linreg(x, y, 1000)

# 计算斜率在95%置信水平下的置信区间
print(np.percentile(bs_slope_reps, [2.5, 97.5]))
[221.48424242 243.84835793]

假设检验

前面我们介绍了参数估计的思路和方法，确定了参数，也就确定了模型，接下来要回答的问题是：如果模型是正确的，数据集能在多大程度上被这个模型所解释？我们可以通过假设检验来回答这个问题。

假设检验的基本思想是小概率反证法。根据所考察问题的要求提出原假设H0和备择假设H1，这两个假设是对立的，其中备择假设H1是我们想要证明正确的结论，整个验证的思路是：在原假设H0成立的条件下，用适当的统计方法证明原假设H0成立的概率很小，进而推出原假设H0不成立，拒绝原假设H0，接受备择假设H1。

比如，原假设为两个变量有完全相同的概率分布。为了检验这个假设，对数据集进行排列抽样是一个很好的方法，我们可以调用np.random.permutation()【3】函数实现对假设的模拟。

模拟的思路是，先将两个数据集合并成一个数组，调用np.random.permutation()函数对该数组进行随机排列，然后再分成两个数据集，如果原假设成立，则随机排列前/后的两组数据集不应该有差别。

df1 = ts.pro_bar(ts_code='000001.SZ', adj='qfq', start_date='20180101', end_date='20180630')
df1.sort_values('trade_date', inplace=True)
df1 = np.array(df1['close'])[:5]
df2 = ts.pro_bar(ts_code='000001.SZ', adj='qfq', start_date='20190101', end_date='20190630')
df2.sort_values('trade_date', inplace=True)
df2 = np.array(df2['close'])[:5]


# 把两个数据集合并成一个数组
data = np.concatenate((df1, df2))
print(data)
[13.7  13.33 13.25 13.3  12.96  9.09  9.18  9.65  9.64  9.56]
# 随机排列抽样
permuted_data = np.random.permutation(data)
print(permuted_data)
[13.25  9.65  9.09  9.56  9.18 13.7  13.33 12.96  9.64 13.3]

# 重新将其分成两个数据集，完成一个模拟。
perm_sample_1 = permuted_data[:len(df1)]
perm_sample_2 = permuted_data[len(df1):]
print(perm_sample_1)
print(perm_sample_2)


[13.7  13.25  9.64 13.33  9.09]
[ 9.56  9.65 13.3  12.96  9.18]

# 封装成一个函数，方便后面进行重复调用。

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

现在有了观测数据和模拟数据，下一步就是要选择合适的检验统计量来计算概率。检验统计量是指在原假设成立的条件下，可以从观测数据和模拟数据集中计算出的一个数值。

这里我们选择两组数据均值的差作为检验统计量，然后进行10000次模拟排列抽样，计算模拟均值差不小于观测均值差的(小概率事件)的概率。

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 diff_of_means(data_1, data_2):
    """计算两个数组的均值差"""

    diff = np.mean(data_1)-np.mean(data_2)

    return diff

# 计算观测数据的均值差
empirical_diff_means = diff_of_means(df1, df2)

# 计算10,000次模拟数据的均值差
perm_replicates = draw_perm_reps(df1, df2, diff_of_means, size=10000)

# 计算P值: p
p = np.sum(perm_replicates >= empirical_diff_means) / len(perm_replicates)

# 打印结果
print('p-value =', p)


5.437731092436975

[-0.40647059 -0.71991597  0.71571429 ...  0.01218487 -0.25773109
  0.40546218]

p-value = 0.0

我们发现，p值为0，说明这个原假设成立的概率很小，所以拒绝H0，接受H1。

假设检验的案例

A/B测试

A/B测试是公司用来检验新策略是否更有效的一种方法。比如，一家公司Web界面有A/B两个版本，其中A是原始版，B是更新版，现在想知道B版本的应用会不会带来更多的点击量，这时就可以用A/B测试。

具体来说，可以将1000名访客分成各500人的两组，随机访问这两个版本的界面，如果结果显示B版本的点击量是67次，A版本是45次，则说明新版本B比A更有效，可以考虑使用。

下面，我们看一个实际的情景。

1920年，美国职业棒球大联盟实施了重要的规则更改，结束了所谓的死球时代，即投手不再被允许往球上吐口水或用手擦球，这是一项非常有利于投手的活动。数据集nht_dead和nht_live分别代表在这项规则使用前/后，两次击中球之间的时间间隔。

我们想要研究的问题是，在这些规则使用之后，是否会导致两次击中球之间的时间间隔变长，所以要进行假设检验。

H0：规则使用后，没有导致两次击中球之间的时间间隔变长。

H1：规则使用后，导致两次击中球之间的时间间隔变长。

同样，我们选择均值差作为检验统计量，具体的模拟过程如下：

nht_dead = np.array([-1,  894,   10,  130,    1,  934,   29,    6,  485,  254,  372,
         81,  191,  355,  180,  286,   47,  269,  361,  173,  246,  492,
        462, 1319,   58,  297,   31, 2970,  640,  237,  434,  570,   77,
        271,  563, 3365,   89,    0,  379,  221,  479,  367,  628,  843,
       1613, 1101,  215,  684,  814,  278,  324,  161,  219,  545,  715,
        966,  624,   29,  450,  107,   20,   91, 1325,  124, 1468,  104,
       1309,  429,   62, 1878, 1104,  123,  251,   93,  188,  983,  166,
         96,  702,   23,  524,   26,  299,   59,   39,   12,    2,  308,
       1114,  813,  887])

nht_live = np.array([645, 2088,   42, 2090,   11,  886, 1665, 1084, 2900, 2432,  750,
       4021, 1070, 1765, 1322,   26,  548, 1525,   77, 2181, 2752,  127,
       2147,  211,   41, 1575,  151,  479,  697,  557, 2267,  542,  392,
         73,  603,  233,  255,  528,  397, 1529, 1023, 1194,  462,  583,
         37,  943,  996,  480, 1497,  717,  224,  219, 1531,  498,   44,
        288,  267,  600,   52,  269, 1086,  386,  176, 2199,  216,   54,
        675, 1243,  463,  650,  171,  327,  110,  774,  509,    8,  197,
        136,   12, 1124,   64,  380,  811,  232,  192,  731,  715,  226,
        605,  539, 1491,  323,  240,  179,  702,  156,   82, 1397,  354,
        778,  603, 1001,  385,  986,  203,  149,  576,  445,  180, 1403,
        252,  675, 1351, 2983, 1568,   45,  899, 3260, 1025,   31,  100,
       2055, 4043,   79,  238, 3931, 2351,  595,  110,  215,    0,  563,
        206,  660,  242,  577,  179,  157,  192,  192, 1848,  792, 1693,
         55,  388,  225, 1134, 1172, 1555,   31, 1582, 1044,  378, 1687,
       2915,  280,  765, 2819,  511, 1521,  745, 2491,  580, 2072, 6450,
        578,  745, 1075, 1103, 1549, 1520,  138, 1202,  296,  277,  351,
        391,  950,  459,   62, 1056, 1128,  139,  420,   87,   71,  814,
        603, 1349,  162, 1027,  783,  326,  101,  876,  381,  905,  156,
        419,  239,  119,  129,  467])


# 计算观测数据集的均值差
nht_diff_obs = diff_of_means(nht_dead, nht_live)

# 进行10,000排列抽样
perm_replicates = draw_perm_reps(nht_dead, nht_live, diff_of_means, size=10000)

# 计算P值
p = np.sum(perm_replicates <= nht_diff_obs) / len(perm_replicates)
print('p-val =', p)

p-val = 0.0002