Skip to content

Commit 672b37e

Browse files
acp29acp29
committed
Added option for bias-corrected intervals to bootint - updated matlab toolbox
1 parent 5872ed7 commit 672b37e

5 files changed

Lines changed: 131 additions & 32 deletions

File tree

‎inst/bootint.m‎

Lines changed: 95 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -1,16 +1,17 @@
1-
% Computes simple percentile confidence interval(s) directly from a vector (or
2-
% row-major matrix) of bootstrap statistics.
1+
% Computes percentile confidence interval(s) directly from a vector (or row-
2+
% major matrix) of bootstrap statistics.
33
%
44
% -- Function File: CI = bootint (BOOTSTAT)
55
% -- Function File: CI = bootint (BOOTSTAT, PROB)
6+
% -- Function File: CI = bootint (BOOTSTAT, PROB, ORIGINAL)
67
%
7-
% 'CI = bootint (BOOTSTAT)' computes 95% percentile confidence intervals
8-
% directly from the vector, or rows* of the matrix in BOOTSTAT, where
9-
% BOOTSTAT contains bootstrap statistics such as those generated using
10-
% the `bootstrp` function. Depending on the application, bootstrap
11-
% confidence intervals with better coverage and accuracy can be computed
12-
% using the various dedicated bootstrap functions from the statistics-
13-
% resampling package.
8+
% 'CI = bootint (BOOTSTAT)' computes simple 95% percentile confidence
9+
% intervals [1,2] directly from the vector, or rows* of the matrix in
10+
% BOOTSTAT, where BOOTSTAT contains bootstrap statistics such as those
11+
% generated using the `bootstrp` function. Depending on the application,
12+
% bootstrap confidence intervals with better coverage and accuracy can
13+
% be computed using the various dedicated bootstrap confidence interval
14+
% functions from the statistics-resampling package.
1415
%
1516
% * The matrix should have dimensions P * NBOOT, where P corresponds to
1617
% the number of parameter estimates and NBOOT corresponds to the number
@@ -28,7 +29,22 @@
2829
% The default value of PROB is the vector: [0.025, 0.975], for an
2930
% equal-tailed 95% percentile confidence interval.
3031
%
31-
% bootint (version 2024.04.24)
32+
% 'CI = bootint (BOOTSTAT, PROB, ORIGINAL)' uses the ORIGINAL estimates
33+
% associated with BOOTSTAT to correct PROB and the resulting confidence
34+
% intervals (CI) for median bias. The confidence intervals returned in CI
35+
% therefore become bias-corrected percentile intervals [3,4].
36+
%
37+
% BIBLIOGRAPHY:
38+
% [1] Efron (1979) Bootstrap Methods: Another look at the jackknife.
39+
% Annals Stat. 7,1-26
40+
% [2] Efron, and Tibshirani (1993) An Introduction to the Bootstrap.
41+
% New York, NY: Chapman & Hall
42+
% [3] Efron (1981) Nonparametric Standard Errors and Confidence Intervals.
43+
% Can J Stat. 9(2):139-172
44+
% [4] Efron (1982) The jackknife, the bootstrap, and other resampling plans.
45+
% SIAM-NSF, CBMS #38
46+
%
47+
% bootint (version 2024.05.19)
3248
% Author: Andrew Charles Penn
3349
% https://www.researchgate.net/profile/Andrew_Penn/
3450
%
@@ -46,13 +62,13 @@
4662
% You should have received a copy of the GNU General Public License
4763
% along with this program. If not, see http://www.gnu.org/licenses/
4864

49-
function CI = bootint (Y, PROB)
65+
function CI = bootint (Y, PROB, T0)
5066

5167
% Check input and output arguments
5268
if (nargin < 2)
5369
PROB = 0.95;
5470
end
55-
if (nargin > 2)
71+
if (nargin > 3)
5672
error ('bootint: Too many input arguments.')
5773
end
5874
if (nargout > 1)
@@ -67,6 +83,7 @@
6783
if (sz(2) == 1)
6884
p = 1;
6985
Y = Y.';
86+
nboot = numel (Y);
7087
else
7188
p = sz(1);
7289
nboot = sz(2);
@@ -100,6 +117,43 @@
100117
PROB = arrayfun (@(c) (1 + c * PROB) / 2, [-1, 1]);
101118
end
102119
end
120+
121+
% Apply median bias correction to PROB if T0 is provided
122+
if (nargin > 2)
123+
% Error checking for T0
124+
if (~ isa (PROB, 'numeric'))
125+
error ('bootint: ORIGINAL must be numeric')
126+
end
127+
szT0 = size (T0);
128+
if ( (all (szT0 > 1)) || (szT0(2) > 1) )
129+
error ('bootint: ORIGINAL must be a scalar or column vector')
130+
end
131+
if (p > 1)
132+
if (szT0(1) ~= p)
133+
error ('bootint: BOOTSTAT and ORIGINAL must have the same number of rows')
134+
end
135+
else
136+
if (szT0(1) > 1)
137+
error ('bootint: If BOOTSTAT is a vector, ORIGINAL must be a scalar')
138+
end
139+
end
140+
% Create distribution functions
141+
stdnormcdf = @(x) 0.5 * (1 + erf (x / sqrt (2)));
142+
stdnorminv = @(p) sqrt (2) * erfinv (2 * p - 1);
143+
% Calculate the median bias correction constant (z0)
144+
z0 = stdnorminv (sum (bsxfun (@lt, Y, T0), 2) / nboot);
145+
if (~ all (isfinite (z0)))
146+
% Revert to percentile bootstrap confidence intervals
147+
warning ('bootint:biasfail', ...
148+
cat (2, 'Unable to calculate the bias correction constant;', ...
149+
' reverting to simple percentile intervals.\n'))
150+
z0 = zeros (p, 1);
151+
end
152+
% Calculate BC percentiles
153+
z = stdnorminv (PROB);
154+
PROB = stdnormcdf (bsxfun (@plus, 2 * z0, z));
155+
end
156+
103157
% Compute confidence intervals
104158
CI = nan (p, 2);
105159
for j = 1:p
@@ -117,12 +171,36 @@
117171

118172
end
119173

174+
%!demo
175+
%!
176+
%! % Law school data
177+
%! data = [576, 3.39; 635, 3.30; 558, 2.81; 578, 3.03; 666, 3.44; ...
178+
%! 580, 3.07; 555, 3.00; 661, 3.43; 661, 3.36; 605, 3.13; ...
179+
%! 653, 3.12; 575, 2.74; 545, 2.76; 572, 2.88; 594, 2.96];
180+
%! x = data(:, 1);
181+
%! y = data(:, 2);
182+
%! r = cor (x, y);
183+
%!
184+
%! % 95% confidence interval for the mean
185+
%! bootstat = bootstrp (4999, @cor, x, y);
186+
%! CI_per = bootint (bootstat,0.95) % 95% simple percentile interval
187+
%! CI_cper = bootint (bootstat,0.95,r) % 95% bias-corrected percentile interval
188+
%!
189+
%! % Please be patient, the calculations will be completed soon...
190+
120191
%!test
121192
%!
122-
%! % Simulate (log-normal) data
123-
%! randn ('seed', 1);
124-
%! Y = exp (randn (5, 999));
193+
%! % Law school data
194+
%! data = [576, 3.39; 635, 3.30; 558, 2.81; 578, 3.03; 666, 3.44; ...
195+
%! 580, 3.07; 555, 3.00; 661, 3.43; 661, 3.36; 605, 3.13; ...
196+
%! 653, 3.12; 575, 2.74; 545, 2.76; 572, 2.88; 594, 2.96];
197+
%! x = data(:, 1);
198+
%! y = data(:, 2);
199+
%! r = cor (x, y);
125200
%!
126201
%! % 95% confidence interval for the mean
127-
%! CI = bootint (Y,0.95); # 95% percentile interval
128-
%! CI = bootint (Y,[0.025,0.975]); # 95% percentile interval
202+
%! bootstat = bootstrp (4999, @cor, x, y);
203+
%! CI = bootint (bootstat,0.95); % 95% percentile interval
204+
%! CI = bootint (bootstat,[0.025,0.975]); % 95% percentile interval
205+
%! CI = bootint (bootstat,0.95,r); % 95% bias-corrected interval
206+
%! CI = bootint (bootstat,[0.025,0.975],r); % 95% bias-corrected interval

‎inst/credint.m‎

Lines changed: 19 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -145,12 +145,26 @@
145145

146146
end
147147

148+
%!demo
149+
%!
150+
%! % Input univariate dataset
151+
%! heights = [183, 192, 182, 183, 177, 185, 188, 188, 182, 185].';
152+
%!
153+
%! % 95% credible interval for the mean
154+
%! [stats, bootstat] = bootbayes (heights);
155+
%! CI = credint (bootstat,0.95) % 95% shortest probability interval
156+
%! CI = credint (bootstat,[0.025,0.975]) % 95% equal-tailed interval
157+
%!
158+
%! % Please be patient, the calculations will be completed soon...
159+
148160
%!test
149161
%!
150-
%! % Simulate (log-normal) data
151-
%! randn ('seed', 1);
152-
%! Y = exp (randn (5, 999));
162+
%! % Input univariate dataset
163+
%! heights = [183, 192, 182, 183, 177, 185, 188, 188, 182, 185].';
164+
%!
165+
%! % 95% credible interval for the mean
166+
%! [stats, bootstat] = bootbayes (heights);
153167
%!
154168
%! % 95% credible interval for the mean
155-
%! CI = credint (Y,0.95); # 95% shortest probability interval
156-
%! CI = credint (Y,[0.025,0.975]); # 95% equal-tailed interval
169+
%! CI = credint (bootstat,0.95); % 95% shortest probability interval
170+
%! CI = credint (bootstat,[0.025,0.975]); % 95% equal-tailed interval

‎matlab/statistics-resampling.mltbx‎

1.09 KB
Binary file not shown.

‎matlab/statistics-resampling.prj‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,13 @@
11
<deployment-project plugin="plugin.toolbox" plugin-version="1.0">
2-
<configuration build-checksum="3022689601" file="Y:\Documents\GitHub\statistics-resampling\matlab\statistics-resampling.prj" location="Y:\Documents\GitHub\statistics-resampling\matlab" name="statistics-resampling" target="target.toolbox" target-name="Package Toolbox">
2+
<configuration build-checksum="2670653735" file="Y:\Documents\GitHub\statistics-resampling\matlab\statistics-resampling.prj" location="Y:\Documents\GitHub\statistics-resampling\matlab" name="statistics-resampling" target="target.toolbox" target-name="Package Toolbox">
33
<param.appname>statistics-resampling</param.appname>
44
<param.authnamewatermark>Andrew Penn</param.authnamewatermark>
55
<param.email>andy.c.penn@gmail.com</param.email>
66
<param.company>University of Sussex, UK</param.company>
77
<param.summary>Statistical analysis using resampling methods</param.summary>
88
<param.description>The statistics-resampling package is an Octave package and Matlab toolbox that can be used to perform a wide variety of statistics tasks using non-parametric resampling methods. In particular, the functions included can be used to estimate bias, uncertainty (standard errors and confidence intervals), prediction error, and calculate p-values for null hypothesis significance tests. Variations of the resampling methods are included that improve the accuracy of the statistics for small samples and samples with complex dependence structures.</param.description>
99
<param.screenshot>Y:\Documents\GitHub\statistics-resampling\doc\icon.png</param.screenshot>
10-
<param.version>5.5.13</param.version>
10+
<param.version>5.5.14</param.version>
1111
<param.output>${PROJECT_ROOT}\statistics-resampling.mltbx</param.output>
1212
<param.products.name />
1313
<param.products.id />

‎test/test_script.m‎

Lines changed: 15 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -764,16 +764,23 @@
764764
pval7 = randtest2 (X, Y, false, [], @(A, B) log (var (A) ./ var (B)), 1);
765765

766766
% credint:test:1
767-
randn ('seed', 1);
768-
Y = exp (randn (5, 999));
769-
CI = credint (Y,0.95); % 95% shortest probability interval
770-
CI = credint (Y,[0.025,0.975]); % 95% equal-tailed interval
767+
heights = [183, 192, 182, 183, 177, 185, 188, 188, 182, 185].';
768+
[stats, bootstat] = bootbayes (heights);
769+
CI_1 = credint (bootstat,0.95); % 95% shortest probability interval
770+
CI_2 = credint (bootstat,[0.025,0.975]); % 95% equal-tailed interval
771771

772772
% bootint:test:1
773-
randn ('seed', 1);
774-
Y = exp (randn (5, 999));
775-
CI = bootint (Y,0.95); % 95% percentile interval
776-
CI = bootint (Y,[0.025,0.975]); % 95% percentile interval
773+
data = [576, 3.39; 635, 3.30; 558, 2.81; 578, 3.03; 666, 3.44; ...
774+
580, 3.07; 555, 3.00; 661, 3.43; 661, 3.36; 605, 3.13; ...
775+
653, 3.12; 575, 2.74; 545, 2.76; 572, 2.88; 594, 2.96];
776+
x = data(:, 1);
777+
y = data(:, 2);
778+
r = cor (x, y);
779+
bootstat = bootstrp (4999, @cor, x, y);
780+
CI_1 = bootint (bootstat,0.95); % 95% percentile interval
781+
CI_2 = bootint (bootstat,[0.025,0.975]); % 95% percentile interval
782+
CI_3 = bootint (bootstat,0.95,r); % 95% bias-corrected interval
783+
CI_4 = bootint (bootstat,[0.025,0.975],r); % 95% bias-corrected interval
777784

778785
% sampszcalc:test:1
779786
ns = sampszcalc ('t', 0.20, 0.80, 0.05, 2);

0 commit comments

Comments
 (0)