BootStrap(自助法):一种有返还的再抽样统计方法,可以用于总体分布未知或统计量的分布未知时的参数推断。再抽样过程可以不依赖于某些经典统计推断假设,如抽样分布的正态性和中心极限定律,而用来计算检验统计量的P值或置信区间的上、下限。
典故:术语“Bootstrap”来自短语“to pull oneself up by one’s bootstraps” (源自西方神话故事“ The Adventures of Baron Munchausen”,男爵掉到了深湖底,没有工具,所以他想到了拎着鞋带将自己提起来)
- 计算机的引导程序boot也来源于此
- 意义:不靠外界力量,而靠自身提升自己的性能,翻译为自助/自举
许多经典统计方法是建立在对样本观察值所在总体进行多种假设基础上的。例如,如果样本呈正态分布,计算
该t值将服从一个自由度为的t分布。然而,在许多情况下,样本观测值是偏离正态分布的活着观测值的分布根本不知道。因此基于正态分布假设的经典统计推断结果便不可靠甚至可能错误。
基本思路:利用样本数据计算统计量和估计样本分布,而不对模型做任何假设(非参数bootstrap)
如果不知道总体分布,那么对总体分布的最好猜测便是由数据提供的分布。
- 假定观测值便是总体
- 由这一假定的总体有放回的再抽样样本,由原始数据经过再抽样所获得的与原始数据集含量相等的样本称为再抽样样本(resamples)或自助样本(bootstrap samples)
- 由原始数据集计算所得的统计量称为观察统计量(observed statistic);由再抽样样本计算所得的统计量称为自助统计量(bootstrap statistic)。其中,“∷”表示二者间的关系,“<=> ”表示等价于。也就是说,通过对自助统计量的研究,就可以了解有关观察统计量与真值的偏离情况。
假设原始数据集包含个观测值,根据有限的观测值可以推测出产生这些数据的模型并找到最优的参数,但是无法知道模型和参数多大程度上可信。因为观测值x通常包含噪声,同时会有一些影响因素没有被纳入模型的考虑范围。
Bootstrap方法通过创造多个“pseudo dataset”来给出模型及参数的可信程度方法如下:
原始数据集包含N个观测值,其中。
- 从Z中有放回地抽出N个点构成第一数据集
- 从Z中有放回地抽出N个点构成第二数据集
- ...
- 从Z中有放回地抽出N个点构成第m数据集
(显然,Z中的点有些重复出现在中,有些点没有出现在中)
- 在第1数据集上训练出模型的一组参数
- 在第2数据集 上训练出模型的一组参数
- ...
- 在第m数据集上训练出模型的一组参数
在一个bootstrap样本集中不包含某个原始样本的概率为
一个bootstrap样本集包含了大约原始样本集的1-0.368 = 0.632,另外0.368的样本没有包括。
比较不同参数作出的预测的variability(方差,可以理解为波动范围大小),在不同样本上训练出来的模型差异越大可信度越小。Efron将这一置信区间估计方法命名为Bootstrap,其含义是不依靠额外的验证数据集(validation dataset), 仅使用自身的原始数据集,给出自身的可信度。
作为一种再抽样方法,自助法由原始数据集重复抽取样本,然后计算一些统计量的值,并用这些数值对一些未知参数的真值做出推断。那么,就某一统计量而言,到底需要多少自助重复呢?然而,当估计置信区间和进行假设检验时,大致需要2 000个自助样本。考虑到目前计算机效率已不是很大的制约因素,建议自助样本数可以为1 000~100 000不等。
Bootstrap Estimate of Mean
将B个自助样本的均数进行平均便得到自助样本均数,用此代表总体均数。令为第个自助样本集的均数,,为自助样本均数,则
Bootstrap Estimate of Standard Error
使用bootstrap估计的标准差。
注意,这里计算的是一组(B个)均数的标准差,即均数的标准误差。因此,不要根据常规计算标准误差的方法再用去除这一标准差。
Efron(1987)给出一个公式以检验不同样本数的效应,并根据他们的经验指出,若估计偏差和方差,通常B=200即可;当B=50时便能提供一些有用的信息。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 |
clc,clear all; load forearm % Sample with replacement from this. % First get the sample size. n = length(forearm); B = 100;% number of bootstrap replicates % Get the value of the statistic of interest. theta = skewness(forearm); % Use unidrnd to get the indices to the resamples. % Note that each column corresponds to indices % for a bootstrap resample. inds = unidrnd(n,n,B);%产生一组离散均匀随机整数 % Extract these from the data. xboot = forearm(inds); % We can get the skewness for each column using the % MATLAB Statistics Toolbox function skewness. thetab = skewness(xboot); seb = std(thetab); % Now show how to do it with MATLAB Statistics Toolbox % function: bootstrp. Bmat = bootstrp(B,'skewness',forearm); % What we get back are the bootstrap replicates. % Get an estimate of the standard error. sebmat = std(Bmat); |
偏度用来衡量分布的不对称程度或偏斜成都。随机变量的偏度是变量的三阶中心矩除以标准差的三次方(n/ (n-1)(n-2)再乘以(样本与均值的偏差的立方和/样本标准差的立方))。
正态分布的偏度为零,若称分布具有负偏离,也称左偏态,此时数据位于均值右边的比位于左边的多;若,呈分布具有正偏离,也称右偏态。当偏度接近0时,可认为分布是对称的。若知道分布有可能在偏度上偏离正态分布时 ,可用偏离来检测分布的正态性。
Bootstrap Estimate of Bias
Bias is defined as the difference between the expected value of the statistic and the parameter:
上式会产生一个less-biased的估计,但是具有更大的方差或标准差。如果估计得到的偏差相比标准差小的多,那么就不要采用这个修正公式。
More bootstrap samples are needed to estimate the bias, than are required to estimate the standard error.Efron(1987) 推荐
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 |
load forearm % Sample with replacement from this. % First get the sample size. n = length(forearm); B = 100;% number of bootstrap replicates % Get the value of the statistic of interest. theta = skewness(forearm); % Use unidrnd to get the indices to the resamples. % Note that each column corresponds to indices % for a bootstrap resample. inds = unidrnd(n,n,B);%产生一组离散均匀随机整数 % Extract these from the data. xboot = forearm(inds); % We can get the skewness for each column using the % MATLAB Statistics Toolbox function skewness. thetab = skewness(xboot); seb = std(thetab); % Use the same replicates from before. % Evaluate the mean using Equation 6.15. meanb = mean(thetab); % Now estimate the bias using Equation 6.17. biasb = meanb - theta; |
Bootstrap Estimate of Confidence Intervals
1)the standard interval
假设检验中的置信区间
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
mu = 45; sig = 15; n = 100; alpha = 0.05; xbar = 47.2; % Get the 95% confidence interval. % Get the value for z_alpha/2 % norminv(probability,mean,standard_dev) %Probability 正态分布的概率值;Mean 分布的算术平均值;Standard_dev 分布的标准偏差. zlo = norminv(1-alpha/2,0,1);%1.9600 zhi = norminv(alpha/2,0,1);%-1.9600 thetalo = xbar - zlo*sig/sqrt(n); thetaup = xbar - zhi*sig/sqrt(n); |
Bootstrap标准置信区间可由下式求得:
其中 is the standard error for the statistic obtained using the bootstrap [Mooney and Duval, 1993].
上式可用于当为正态分布或正态假设时。
2)the bootstrap-t interval
注意观察上式,为重抽样集的估计标准差。那么,当我们放回重采样一个数据集时可以得到的一个估计。那么我们可以得到个的bootstrap估计。
6.23式说明估计的分位数(quantile)为使得的点将小于该数目。比如,,,那么将被估计为第五大的值。
where is an estimate of the standard error of .
注意,如果无法计算的标准差公式,那么需要采用Bootstrap方法,即意味着有两成bootstrapping:一个用来计算,一个用来计算。比如B=1000,采用50个bootstrap重采样集来寻找算,那么总共需要50000次重采样。
例子:估计forearm数据集的方差即二阶中心矩。
1 2 3 4 5 6 7 8 9 10 11 |
function mr = mom(x) % MOM Sample second central moment. % MR = MOM(X) % This function returns the sample second central moment. n = length(x); mu = mean(x); mr = (1/n)*sum((x-mu).^2); |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 |
clc,clear all; load forearm n = length(forearm); alpha = 0.1; B = 1000; thetahat = mom(forearm); % Get the bootstrap replicates and samples. [bootreps, bootsam] = bootstrp(B,'mom',forearm); % Set up some storage space for the SE抯. sehats = zeros(size(bootreps)); % Each column of bootsam contains indices % to a bootstrap sample. for i = 1:B % extract the sample from the data xstar = forearm(bootsam(:,i)); bvals(i) = mom(xstar); % Do bootstrap using that sample to estimate SE. sehats(i) = std(bootstrp(25,'mom',xstar)); end zvals = (bootreps - thetahat)./sehats; % Estimate the SE using the bootstrap. SE = std(bootreps); % Get the quantiles. k = B*alpha/2; szval = sort(zvals); tlo = szval(k); thi = szval(B-k); % Get the endpoints of the interval. blo = thetahat-thi*SE; bhi = thetahat-tlo*SE; |
3)the percentile method
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
% Example 7.12 % Computational Statistics Handbook with MATLAB, 2nd Edition % Wendy L. and Angel R. Martinez % Use Statistics Toolbox function % to get the bootstrap replicates. bvals = bootstrp(B,'mom',forearm); % Find the upper and lower endpoints k = B*alpha/2; sbval = sort(bvals); blo = sbval(k); bhi = sbval(B-k); |
置信区间为(1.03, 1.45),相对于第二种方法置信区间较窄。
讨论:
- The standard interval is the easiest and assumes that is normally distributed.
- The bootstrap-t interval estimates the standardized version of from the data, avoiding the normality assumptions used in the standard interval.
- The per-centile interval is simple to calculate and obtains the endpoints directly from the bootstrap estimate of the distribution for . It has another advantage in that it is range-preserving. This means that if the parameter can take on values in a certain range, then the confidence interval will reflect that. This is not always the case with the other intervals.
According to Efron and Tibshirani [1993], the bootstrap-t interval has good coverage probabilities, but does not perform well in practice. The bootstrap percentile interval is more dependable in most situations, but does not enjoy the good coverage property of the bootstrap-t interval. There is another boot-strap confidence interval, called the BCa interval, that has both good cover-age and is dependable. This interval is described in the next chapter.
附录
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 |
function [mu, bias, varjack] = mean_bootstrap(data, B) %Find the estimate of the mean, it's bias and variance using the bootstrap estimator method %Inputs: % data - The data from which to estimate % B - Number of sets to draw %Outputs: % mu - The mean % bias - The bias of the estimator % var - The variance of the estimate [D, N] = size(data); mu_star = zeros(D,B); for i = 1:B, %Draw N samples from the data, with replacement indices = zeros(1,N); for j = 1:N, indices(j) = 1 + floor(rand(1)*N); end %Average the data points chosen mu_star(:,i) = mean(data(:,indices)')'; end mu = 1/B*sum(mu_star')'; bias = mu-mean(data')'; varjack = 1/B*sum((mu_star-mu*ones(1,B))'.^2)'; |