bstrapdemo.m

来自「bootstrp 抽样方法」· M 代码 · 共 46 行

% BSTRAPDEMO  Demonstration of BSTRAP and BSTATS
% EXAMPLE   : bstrapdemo (Oh, the FEX code metrics..)  
% NOTES     : Function HISTF needs to be downloaded from FEX
% AUTHOR    : Dimitri Shvorob, [e-mail], 4/15/07

B = 1000;

% Simple regression
n = 100;
x = rand(n,1);
i = ones(n,1);
a = 3;
b = 5;
y = a + b*x + .5*randn(n,1);

% Bootstrap calculations
disp('A. "fun" with vector output (OLS estimates of regression intercept and slope)')
ls = bstrap(B,1,'olsvector',[y i x]); 
disp('Mean and variance estimates, plus two-sided 95% confidence intervals')
[m,v,cil,ciu] = bstats(ls)        %#ok  note: may want to exclude actual-sample value 

disp('B. "fun" with structure output (OLS estimates of regression intercept and slope)')
ls = bstrap(B,1,'olsstructure',[y i x]);   
disp('Mean and variance estimates, plus two-sided 95% confidence intervals')
[m,v,cil,ciu] = bstats(ls,'b')    %#ok  note: may want to exclude actual-sample value 

% Bootstrap distributions of OLS estimates
ab = [ls.b];                      % note: may want to exclude actual-sample value 
subplot(3,1,1)
scatter(x,y,2)
title(['Random sample: \ity = ' num2str(a) ' + ' num2str(b) 'x + \epsilon, \rmwhere \itx ~ U\rm[\it0,1\rm]\it,\epsilon ~ N\rm[\it0,1\rm], \iti = 1,..,100']); 

subplot(3,1,2)
opt.dx   = .01;
opt.xmin = a - .5;
opt.xmax = a + .5;
opt.xmrk = a; 
histf(ab(1,:),opt)
title(['OLS intercept estimate: bootstrap distribution. (True value = \it' num2str(a) '\rm)'])

subplot(3,1,3)
opt.xmrk = b;
opt.xmin = b - .5;
opt.xmax = b + .5;
histf(ab(2,:),opt)
title(['OLS slope estimate: bootstrap distribution. (True value = \it' num2str(b) '\rm)'])