2.5 模型可信度

一、有放回和无放回抽样

当我们进行模型拟合时，模型中的参数通常都需要估计，这个参数估计的结果（通常体现在误差棒上）能直接体现模型是否可靠。不同的抽样方法采用的抽样方法不同，但主要分为有放回和无放回抽样。以在一个有N盒子中抽取f个球为例：

有放回地抽样（sample with replacement）指每次抽取1个球，然后放回去，重复f次，总共有N^f种可能，可能取到相同的数据；
无放回地抽样（sample without replacement）指直接抽取f个球，总共有C(N, f)种可能（在不考虑抽取球顺序的情况下），不会抽到相同的数据。

二、自举法（bootstrap）

自举法（bootstrap）是一种常用的参数区间估计方法，采用有放回地抽样。

举个例子，对于一个有N个样本的数据集，模型的形式已经确定：y=ax+b，(a,b)为需要估计的参数。自举法则从这些数据中有放回地抽取N个数据，根据抽取的数据对模型的参数进行估计，得到一组参数(a,b)；如果重复抽取1000次（bootstrap=1000），则会得到1000组参数(a,b)。根据这些估计的结果，我们可以得到参数a和b的分布。当我们取分布中的2.5%至97.5%范围，就可以得到95%CI（也常用作绘制误差棒error bars），它的物理意义是，我们认为真实的参数a（或b）有95%的可能落在该区间。换句话说，我们有95%的信心认为这个区间包括了真实的参数值。容易看出，这个区间越窄，我们对于参数的估计值就越有自信。

举个例子，在高考后需要估分填志愿的情况下，如果一个人估计自己的高考成绩均值为600分，并认为自己的高考成绩有95%的可能在400分至700分之间，这对志愿填报几乎没有帮助；而同样是均值为600分，如果一个人认为自己的高考成绩有95%的可能在590至610之间，那么这对于志愿填报有相当大的参考意义，几乎可以锁定学校的层次。

因此，我们希望得到一个尽可能窄的参数估计区间。一般情况下，增加数据集的样本量N能缩小参数估计区间。图2.5.1给出了一个估计区间受到样本量影响的例子。上下两份相似的图像中，其真模型的形式均为y=0.5x+2，拟合模型的形式均为y=ax+b，其唯一的区别在于采用自举法时采集的样本量不同，上半部分的样本量为n=10，下半部分为n=200，这就导致参数a与b分布的估计区间宽窄不同（见右边条形图），上方的红色实线为95%CI（误差棒）。

可以看出，当样本量很少时（n=10），对于参数a和b的估计区间较宽，误差棒也较大，均值也和真实模型有较大的差距；而当样本量较大时（n=200），对于参数a和b的估计区间很窄，误差棒也很小，几乎接近真实模型。

通过代码，我们来模拟使用自举法来估计模型可靠性：

我们先假定真实模型 $y=ax+b$ , $a = 0.5$ , $b = 3$ 来生成数据：

import numpy as np
import matplotlib.pyplot as plt 

a = 0.5 # coefficient
b = 3   # intercept
noiseVar = 5 # noise的标准差
x = np.arange(1, 100, 1)
# generate data
y =  a * x + b + noiseVar* np.random.randn(x.size)

# 画出数据的散点图
plt.plot(x, y, 'o', color='r')
plt.xlabel('Population')
plt.ylabel('House price')

可以得到以下分布图：

我们用bootstrap的方法来估计参数的可靠性(reliability)

nBoot = 1000 # bootstrap多少次
res = np.empty((nBoot, 2))
for i in range(nBoot):
    # sample with replacement data
    idx = np.random.choice(np.arange(y.size), y.size, replace=True)

    # deal with x
    x2 = np.vstack((x[idx], np.ones(x.size))) # 注意这里x的idx也要改变

    # fit the square
    res[i, :] = np.linalg.lstsq(x2.T, y[idx], rcond=None)[0]

绘制分布图：

# 画出分布图
fig, ax = plt.subplots(1, 3, figsize=(12, 4))
# subplot 1
plt.sca(ax[0])
plt.plot(x, y, 'o', color='r')
plt.xlabel('Independent variable')
plt.ylabel('Dependent variable')
for i in range(nBoot):
    plt.plot(x, res[i,0]*x+res[i,1], '-', color='gray')

# subplot 2
plt.sca(ax[1])
plt.hist(res[:, 0])
plt.xlim([0.4, 0.6])
plt.xlabel('Coefficient')
plt.ylabel('Count')

# subplot 3
plt.sca(ax[2])
plt.hist(res[:, 1])
plt.xlim([0, 10])
plt.xlabel('Intercept')
plt.ylabel('Count')

得到分布图

如果我们有更多的数据，那么估计精度会出现什么样的变化呢？我们首先生成更多的模拟数据：

import numpy as np
import matplotlib.pyplot as plt 

a = 0.5 # coefficient
b = 3   # intercept
noiseVar = 5 # noise的标准差
x = np.arange(1, 1000, 0.1)
# generate data
y =  a * x + b + noiseVar* np.random.randn(x.size)
plt.plot(x, y, 'o', color='r')

得到的散点图为：

现在用原来10倍的数据，再来用bootstrap来估计原来的reliability

nBoot = 1000 # bootstrap多少次
res = np.empty((nBoot, 2))
for i in range(nBoot):
    # sample with replacement data
    idx = np.random.choice(np.arange(y.size), y.size, replace=True)

    # deal with x
    x2 = np.vstack((x[idx], np.ones(x.size))) # 注意这里x的idx也要改变

    # fit the square
    res[i, :] = np.linalg.lstsq(x2.T, y[idx], rcond=None)[0]

# 画出分布图
fig, ax = plt.subplots(1, 3, figsize=(12, 4))
# subplot 1
plt.sca(ax[0])
plt.plot(x, y, 'o', color='r')
plt.xlabel('Independent variable')
plt.ylabel('Dependent variable')
for i in range(nBoot):
    plt.plot(x, res[i,0]*x+res[i,1], '-', color='gray')

# subplot 2
plt.sca(ax[1])
plt.hist(res[:, 0])
plt.xlim([0.4, 0.6])
plt.xlabel('Coefficient')
plt.ylabel('Count')

# subplot 3
plt.sca(ax[2])
plt.hist(res[:, 1])
plt.xlim([0, 10])
plt.xlabel('Intercept')
plt.ylabel('Count')

得到新的分布图

综合两次自举法的结果，可以发现，增加数据集的样本量N能缩小参数估计区间，也就是得到更可靠的参数估计区间。

三、模型可靠性和准确性的关系

模型可靠性和准确性都是评价模型好坏的重要指标，可靠性与准确性的关系和信度与效度的关系相似。正如效度建立在信度的基础之上，模型准确性也建立在可靠性的基础上。如图2.3.7所示，黑色曲线为真实模型，蓝色曲线为拟合模型，浅蓝色的区域为拟合模型的95%CI（或误差棒）。当可靠性很高（浅蓝色区域很窄）时，模型准确性可能很高（左上），也可能很低（右上）；当可靠性很低（浅蓝色区域很宽）时，模型准确性自然也很低了（右下）。一般不会出现可靠性很低但准确性很高的情况。

理论上，自举法和交叉验证几乎使用于所有模型，即通过自举法评估模型可靠性，通过交叉验证评估模型准确性（如图2.5.7所示）。考虑到计算机算力有限，研究者可以根据实际情况选取合适的自举法次数和交叉验证的具体过程。

四、置换检验

建模中另一个重要的问题是，如何知道我们估计出来的参数是显著的？置换检验(premutation test) 是一种常用的检验方法，它遵循假设检验的逻辑。举个例子，我们想要拟合人口(x)和房价(y)的关系，确定了模型的形式为y=ax+b。那么，零假设为人口和房价没有任何关系，即a=0；备择假设为人口和房价有关系，即a≠0。假如我们得到a=0.45，如何确定它显著大于0呢？我们可以使用置换检验的方法，将数据点随机打乱（即无放回地抽样），重新估计参数a的值，并重复多次，则可以得到随机数据后的参数a的分布。将0.45和该分布进行比较，如果处于显著的位置（例如0.05的显著性水平），则可以认为a=0.45是显著大于0的。

置换检验是一种思想，可以用到很多的场景和参数检验。基本流程是构建一个某个参数的零假设的分布，然后利用得到的参数值和该分布做对比。如果得到的参数值在这个分布中处于显著性水平（例如0.05）的两端，则表示该参数值在随机情况下出现的概率较低，从而可以认为该参数值是显著的。

利用代码模拟置换检验，以得到线性系数的显著性p-value：

我们先假定真实模型 $y=ax+b$ , $a = 0.5$ , $b = 3$ 来生成数据：

import numpy as np
import matplotlib.pyplot as plt 

a = 0.05 # coefficient
b = 3   # intercept
noiseVar = 5 # noise的标准差
x = np.arange(1, 100, 0.1)
# generate data
y =  a * x + b + noiseVar* np.random.randn(x.size)

然后我们来fit一条直线

x2 = np.vstack((x, np.ones(x.size))) # 注意这里的x不用利用idx而改变
res = np.linalg.lstsq(x2.T, y, rcond=None)[0]
a = res[0]
b = res[1]
print(res)

得到（a,b）的结果是：

[0.03829224 3.6987213 ]

我们用permutation的方法来得到线性系数的显著性p-value

nBoot = 10000 # permutation 多少次
res2 = np.empty((nBoot, 2))
for i in range(nBoot):
    # sample without replacement data
    idx = np.random.choice(np.arange(y.size), y.size, replace=False)

    # deal with x
    x2 = np.vstack((x, np.ones(x.size))) # 注意这里的x不用利用idx而改变

    # fit the square
    res2[i, :] = np.linalg.lstsq(x2.T, y[idx], rcond=None)[0]

绘制分布图：

# 画出分布图
plt.hist(res2[:, 0])
plt.xlabel('Coefficient')
plt.ylabel('Count')

得到参数a的零假设分布图为

计算并输出p-value

print('Significance probability is ', (res2[:, 0]>a).sum()/nBoot)

结果为

Significance probability is  0.0

可以看出，p<0.05，参数a是显著的。

最后更新于5个月前

import numpy as np import matplotlib.pyplot as plt a = 0.5 # coefficient b = 3 # intercept noiseVar = 5 # noise的标准差 x = np.arange(1, 100, 1) # generate data y = a * x + b + noiseVar* np.random.randn(x.size) # 画出数据的散点图 plt.plot(x, y, 'o', color='r') plt.xlabel('Population') plt.ylabel('House price')

nBoot = 1000 # bootstrap多少次 res = np.empty((nBoot, 2)) for i in range(nBoot): # sample with replacement data idx = np.random.choice(np.arange(y.size), y.size, replace=True) # deal with x x2 = np.vstack((x[idx], np.ones(x.size))) # 注意这里x的idx也要改变 # fit the square res[i, :] = np.linalg.lstsq(x2.T, y[idx], rcond=None)[0]

# 画出分布图 fig, ax = plt.subplots(1, 3, figsize=(12, 4)) # subplot 1 plt.sca(ax[0]) plt.plot(x, y, 'o', color='r') plt.xlabel('Independent variable') plt.ylabel('Dependent variable') for i in range(nBoot): plt.plot(x, res[i,0]*x+res[i,1], '-', color='gray') # subplot 2 plt.sca(ax[1]) plt.hist(res[:, 0]) plt.xlim([0.4, 0.6]) plt.xlabel('Coefficient') plt.ylabel('Count') # subplot 3 plt.sca(ax[2]) plt.hist(res[:, 1]) plt.xlim([0, 10]) plt.xlabel('Intercept') plt.ylabel('Count')

import numpy as np import matplotlib.pyplot as plt a = 0.5 # coefficient b = 3 # intercept noiseVar = 5 # noise的标准差 x = np.arange(1, 1000, 0.1) # generate data y = a * x + b + noiseVar* np.random.randn(x.size) plt.plot(x, y, 'o', color='r')

import numpy as np import matplotlib.pyplot as plt a = 0.05 # coefficient b = 3 # intercept noiseVar = 5 # noise的标准差 x = np.arange(1, 100, 0.1) # generate data y = a * x + b + noiseVar* np.random.randn(x.size)

nBoot = 10000 # permutation 多少次 res2 = np.empty((nBoot, 2)) for i in range(nBoot): # sample without replacement data idx = np.random.choice(np.arange(y.size), y.size, replace=False) # deal with x x2 = np.vstack((x, np.ones(x.size))) # 注意这里的x不用利用idx而改变 # fit the square res2[i, :] = np.linalg.lstsq(x2.T, y[idx], rcond=None)[0]