-
Notifications
You must be signed in to change notification settings - Fork 0
/
ex7_ch3.m
95 lines (92 loc) · 3.35 KB
/
ex7_ch3.m
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
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
% Chrysa Tsimperi
% Data Analysis 2021
% Chapter 3 Excerise 7
% Comparison of percentile bootstrap and parametric confidence interval of
% mean.
% Note that the Matlab default bootstrap (bias corrected and accelerated
% percentile method) is also called in for comparison.
clear all
n = 10;
mu = 0;
sigma = 1;
M = 100;
B = 1000;
alpha = 0.05;
dosquare = 0; % If 1, do square transform first
nbin = 10;
% rng(1);
% Generation of all M samples
xM = mu*ones(n,M) + sigma*ones(n,M).*randn(n,M);
if dosquare
xM = xM.^2;
mu = 1; % This is the mean of the Chi-square distribution with one degree of freedom.
end
xmV = mean(xM)';
xsdV = std(xM)';
cimxM = NaN(M,6); % the two first are the limits of parametric ci for mean
% the other two for the percentile bootstrap ci for mean,
% the final two for the bias corrected and accelerated
% percentile method, which is the default bootstrap in Matlab
klower = floor((B+1)*alpha/2);
kup = B+1-klower;
tailpercV = [klower kup]*100/B;
tcrit = tinv(1-alpha/2,n-1);
for iM = 1:M
% Parametric confidence interval (ci) for mean
mx1err = tcrit * xsdV(iM) / sqrt(n);
cimxM(iM,1:2) = [xmV(iM)-mx1err xmV(iM)+mx1err];
% Percentile bootstrap ci for mean
bootmxV = NaN(B,1);
for iB=1:B
rV = unidrnd(n,n,1);
xbV = xM(rV,iM);
bootmxV(iB) = mean(xbV);
end
cimxM(iM,3:4) = prctile(bootmxV,tailpercV);
cimxM(iM,5:6) = bootci(B,@mean,xM(:,iM));
end
% For min
mxmin = min(min(cimxM(:,[1 3 5])));
mxmax = max(max(cimxM(:,[1 3 5])));
xcenterV = linspace(mxmin,mxmax,nbin+1);
xcenterV = xcenterV(1:nbin)+0.5*(xcenterV(2)-xcenterV(1));
parNx = hist(cimxM(:,1),xcenterV);
bootNx = hist(cimxM(:,3),xcenterV);
boot2Nx = hist(cimxM(:,5),xcenterV);
figure(1)
clf
plot(xcenterV,parNx/M,'.-','linewidth',1.5)
hold on
plot(xcenterV,bootNx/M,'.-r','linewidth',1.5)
plot(xcenterV,boot2Nx/M,'.-g','linewidth',1.5)
xlabel('lower ci limit')
ylabel('relative frequency')
title(sprintf('M=%d, B=%d n=%d, lower limit of ci of mean',M,B,n))
legend('parametric','Perc.bootstrap','BiasCorr.bootstrap','Location','Best')
% For max
mxmin = min(min(cimxM(:,[2 4 6])));
mxmax = max(max(cimxM(:,[2 4 6])));
xcenterV = linspace(mxmin,mxmax,nbin+1);
xcenterV = xcenterV(1:nbin)+0.5*(xcenterV(2)-xcenterV(1));
parNx = hist(cimxM(:,2),xcenterV);
bootNx = hist(cimxM(:,4),xcenterV);
boot2Nx = hist(cimxM(:,6),xcenterV);
figure(2)
clf
plot(xcenterV,parNx/M,'.-','linewidth',1.5)
hold on
plot(xcenterV,bootNx/M,'.-r','linewidth',1.5)
plot(xcenterV,boot2Nx/M,'.-g','linewidth',1.5)
xlabel('upper ci limit')
ylabel('relative frequency')
title(sprintf('M=%d, B=%d n=%d, upper limit of ci of mean',M,B,n))
legend('parametric','Perc.bootstrap','BiasCorr.bootstrap','Location','Best')
iparV = find(cimxM(:,1)<mu & cimxM(:,2)>mu);
ibootV = find(cimxM(:,3)<mu & cimxM(:,4)>mu);
iboot2V = find(cimxM(:,5)<mu & cimxM(:,6)>mu);
fprintf('Coverage probability of %1.2f%% parametric ci of mu=%1.2f= %1.3f \n',...
(1-alpha)*100,mu,length(iparV)/M);
fprintf('Coverage probability of %1.2f%% percentile bootstrap ci of mu=%1.2f= %1.3f \n',...
(1-alpha)*100,mu,length(ibootV)/M);
fprintf('Coverage probability of %1.2f%% bias corrected and accelerated percentile bootstrap ci of mu=%1.2f= %1.3f \n',...
(1-alpha)*100,mu,length(iboot2V)/M);