天天看點

自助法 matlab,Bootstrap 自助法

Bootstrap 自助法

2018-12-27

一、簡介部落格欄目

Bootstrap是一個很通用的算法,用來估計标準誤差、置信區間和偏差。由Bradley Efron于1979年提出,用于計算任意估計的标準誤差。術語“Bootstrap”來自短語“to pull oneself up by one’s bootstraps” (源自西方神話故事“TheAdventures of Baron Munchausen”,男爵掉到了深湖底,沒有工具,是以他想到了拎着鞋帶将自己提起來)。

在統計學中,自助法(BootstrapMethod,Bootstrapping或自助抽樣法)可以指任何一種有放回的均勻抽樣,也就是說,每當選中一個樣本,它等可能地被再次選中并被再次添加到訓練集中。自助法能對采樣估計的準确性(标準誤差、置信區間和偏差)進行比較好的估計,它基本上能夠對任何采樣分布的統計量進行估計。

Bootstrap有兩種形式:非參數bootstrap和參數化的bootstrap,但基本思想都是模拟。參數化的bootstrap假設總體的分布已知或總體的分布形式已知,可以由樣本估計出分布參數,再從參數化的分布中進行再采樣,類似于MC。非參數化的bootstrap是從樣本中再抽樣,而不是從分布函數中進行再抽樣。

二、非參數化Bootstrap

基本思想是:假設是我們的估計量為

自助法 matlab,Bootstrap 自助法

,樣本大小為N,從樣本中有放回的再抽樣N個樣本,原來每一個樣本被抽中的機率相同,均為1/N,得到新的樣本我們稱為Bootstrap樣本,重複B次之後我們得到B個bootstrap樣本集,在每一個樣本集上都有對應的估計量

自助法 matlab,Bootstrap 自助法

,對于B個

自助法 matlab,Bootstrap 自助法

,我們可以計算得到标準誤,置信區間,偏置等。

自助法 matlab,Bootstrap 自助法

三、參數化Bootstrap

和非參數化Bootstrap不同的地方在于總體分布函數的形式是已知的,需要根據樣本估計參數  (參數估計),這樣得到經驗分布函數,從經驗分布函數中再采樣得到Bootstrap樣本,非參數化Bootstrap是從原始樣本中再抽樣,得到的Bootstrap樣本與原始樣本有重合。

四、Matlab執行個體

假設我們的整體(population)來自與Bernouli分布(擲硬币),參數theta等于0.7,即投一次有0.7的機率出現1。為了考察采樣點對估計的影響,我們分别采樣了10和100個樣本,采用了參數和非參數方法。

[plain]

自助法 matlab,Bootstrap 自助法
自助法 matlab,Bootstrap 自助法

%%  Bootstrap demo for the MLE for a Bernoulli

close all;

clear all;

%統計量或估計量,這裡是均值

estimator = @mean;

%Bernoulli分布的參數theta

theta = 0.7;

%樣本數目分别為10和100

Ns = [10 100];

for Ni=1:length(Ns)

%目前的樣本數

N = Ns(Ni);

%再采樣次數B

B = 2000;

%參數為theta的Bernoulli分布的采樣

X = rand(1,N) 

%MLE參數估計

bmle = estimator(X);

%參數化Bootstrap方法對應的B次再采樣的統計量

mleBoot = zeros(1,B);

%非參數化Bootstrap方法對應的B次再采樣的統計量

mleBootNP=zeros(1,B);

for b=1:B

%參數化Bootstrap對應的分布函數就是Bernoulli分布,其參數為Mle估計的參數,參數化Bootstrap再采樣得到樣本

Xb = rand(1,N) 

%對再采樣的Bootstrap樣本求統計量

mleBoot(b) = estimator(Xb);

%非參數化Bootstrap方法再采樣,從原來的樣本中再采樣

ndx = unidrnd(N,1,N);

Xnonparam = X(ndx);

%求統計量

mleBootNP(b) = estimator(Xnonparam);

end

%% 繪圖

%參數化Bootstrap繪圖

figure;

hist(mleBoot)

set(gca,"xlim",[0 1]);

%标準誤

se=std(mleBoot)/sqrt(B);

ttl = sprintf("Boot: true = %3.2f, n=%d, mle = %3.2f, se = %5.3f\n", ...

theta, N,mean(mleBoot), se);

title(ttl);

%非參數化Bootstrap繪圖

figure;

hist(mleBootNP)

set(gca,"xlim",[0 1])

nonParaSe=std(mleBootNP)/sqrt(B);

ttl = sprintf("NP boot: true = %3.2f, n=%d, mle = %3.2f, se = %5.3f\n", theta, N, mean(mleBootNP),nonParaSe);

title(ttl);

end

[plain]

function [out] = bootstrap(data,B)

%

%   Bootstrap method -produce B dataSet

%

%   Inputs:

%       data:origianl data set whose size is [M,N],which means feature Dimens is M and the original dataset contains N samples

%       B:number of dataSet

%    Outputs:

%       out:B bootstrap resamples of the input data  whose size is [M,N,B]

[M,N]=size(data);

% by default B=N;

if (exist("B")~=1), B=N;  end;

out=zeros(M,N,B);

index=unidrnd(N,N,B);

out=reshape(data(:,index),M,N,B);

end

對于參數化Bootstrap方法,假設我們已知總體是Bernouli分布,先從樣本中做MLE估計,得出參數theta,這樣我們就可以從分布函數中直接抽樣,而不是像非參數Bootstrap一樣從樣本中再采樣。

to be continued....

免責聲明:本文僅代表文章作者的個人觀點,與本站無關。其原創性、真實性以及文中陳述文字和内容未經本站證明,對本文以及其中全部或者部分内容文字的真實性、完整性和原創性本站不作任何保證或承諾,請讀者僅作參考,并自行核實相關内容。

http://www.pinlue.com/style/images/nopic.gif