Uses bootstrap to evaluate the likely number of real peaks (i.e. modes)
in the distribution of a single set of data.
-- Function File: H = bootmode (X, M)
-- Function File: H = bootmode (X, M, NBOOT)
-- Function File: H = bootmode (X, M, NBOOT, KERNEL)
-- Function File: H = bootmode (X, M, NBOOT, KERNEL, NPROC)
-- Function File: [H, P] = bootmode (X, M, ...)
-- Function File: [H, P, CRITVAL] = bootmode (X, M, ...)
'H = bootmode (X, M)' tests whether the distribution underlying the
univariate data in vector X has M modes. The method employs the
smooth bootstrap as described [1]. The parsimonious approach is to
iteratively call this function, each time incrementally increasing
the number of modes until the null hypothesis (H0) is accepted (i.e.
H=0), where H0 corresponds to the number of modes being equal to M.
- If H = 0, H0 cannot be rejected at the 5% significance level.
- If H = 1, H0 can be rejected at the 5% significance level.
'H = bootmode (X, M, NBOOT)' sets the number of bootstrap replicates
'H = bootmode (X, M, NBOOT, KERNEL)' sets the kernel for kernel
density estimation. Possible values are:
o 'Gaussian' (default)
o 'Epanechnikov'
'H = bootmode (X, M, NBOOT, KERNEL, NPROC)' sets the number of parallel
processes to use to accelerate computations. This feature requires the
Parallel package (in Octave), or the Parallel Computing Toolbox (in
Matlab).
'[H, P] = bootmode (X, M, ...)' also returns the two-tailed p-value of
the bootstrap hypothesis test.
'[H, P, CRITVAL] = bootmode (X, M, ...)' also returns the critical
bandwidth (i.e.the smallest bandwidth achievable to obtain a kernel
density estimate with M modes)
Bibliography:
[1] Efron and Tibshirani. Chapter 16 Hypothesis testing with the
bootstrap in An introduction to the bootstrap (CRC Press, 1994)
bootmode (version 2023.05.02)
Author: Andrew Charles Penn
https://www.researchgate.net/profile/Andrew_Penn/
Copyright 2019 Andrew Charles Penn
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see http://www.gnu.org/licenses/
The following code
% Stamp data example used in reference [1] in bootstrap R package
x=[0.060;0.064;0.064;0.065;0.066;0.068;0.069;0.069;0.069;0.069;0.069; ...
0.069;0.069;0.070;0.070;0.070;0.070;0.070;0.070;0.070;0.070;0.070; ...
0.070;0.070;0.070;0.070;0.070;0.070;0.070;0.070;0.070;0.070;0.070; ...
0.070;0.070;0.070;0.070;0.070;0.070;0.071;0.071;0.071;0.071;0.071; ...
0.071;0.071;0.071;0.071;0.071;0.071;0.071;0.071;0.071;0.071;0.071; ...
0.071;0.071;0.071;0.071;0.072;0.072;0.072;0.072;0.072;0.072;0.072; ...
0.072;0.072;0.072;0.072;0.072;0.072;0.072;0.072;0.072;0.072;0.072; ...
0.072;0.072;0.072;0.072;0.072;0.072;0.072;0.072;0.072;0.072;0.072; ...
0.072;0.072;0.072;0.073;0.073;0.073;0.073;0.073;0.073;0.073;0.073; ...
0.073;0.073;0.073;0.074;0.074;0.074;0.074;0.074;0.074;0.074;0.074; ...
0.074;0.074;0.075;0.075;0.075;0.075;0.075;0.075;0.075;0.075;0.075; ...
0.075;0.075;0.075;0.075;0.075;0.075;0.075;0.075;0.075;0.075;0.075; ...
0.076;0.076;0.076;0.076;0.076;0.076;0.076;0.076;0.076;0.076;0.076; ...
0.076;0.076;0.076;0.076;0.076;0.076;0.076;0.077;0.077;0.077;0.077; ...
0.077;0.077;0.077;0.077;0.077;0.077;0.077;0.078;0.078;0.078;0.078; ...
0.078;0.078;0.078;0.078;0.078;0.078;0.078;0.078;0.078;0.078;0.078; ...
0.078;0.078;0.078;0.078;0.078;0.078;0.078;0.078;0.079;0.079;0.079; ...
0.079;0.079;0.079;0.079;0.079;0.079;0.079;0.079;0.079;0.079;0.079; ...
0.079;0.079;0.079;0.079;0.079;0.079;0.079;0.079;0.079;0.079;0.079; ...
0.079;0.079;0.079;0.079;0.079;0.079;0.079;0.079;0.079;0.079;0.079; ...
0.079;0.079;0.079;0.079;0.079;0.079;0.080;0.080;0.080;0.080;0.080; ...
0.080;0.080;0.080;0.080;0.080;0.080;0.080;0.080;0.080;0.080;0.080; ...
0.080;0.080;0.080;0.080;0.080;0.080;0.080;0.080;0.080;0.080;0.080; ...
0.080;0.080;0.080;0.080;0.080;0.080;0.080;0.080;0.080;0.080;0.081; ...
0.081;0.081;0.081;0.081;0.081;0.081;0.081;0.081;0.081;0.081;0.081; ...
0.081;0.081;0.081;0.082;0.082;0.082;0.082;0.082;0.082;0.082;0.082; ...
0.082;0.082;0.082;0.082;0.082;0.082;0.082;0.082;0.082;0.082;0.083; ...
0.083;0.083;0.083;0.083;0.083;0.083;0.084;0.084;0.084;0.085;0.085; ...
0.086;0.086;0.087;0.088;0.088;0.089;0.089;0.089;0.089;0.089;0.089; ...
0.089;0.089;0.089;0.089;0.090;0.090;0.090;0.090;0.090;0.090;0.090; ...
0.090;0.090;0.091;0.091;0.091;0.092;0.092;0.092;0.092;0.092;0.093; ...
0.093;0.093;0.093;0.093;0.093;0.094;0.094;0.094;0.095;0.095;0.096; ...
0.096;0.096;0.097;0.097;0.097;0.097;0.097;0.097;0.097;0.098;0.098; ...
0.098;0.098;0.098;0.099;0.099;0.099;0.099;0.099;0.100;0.100;0.100; ...
0.100;0.100;0.100;0.100;0.100;0.100;0.100;0.100;0.100;0.100;0.100; ...
0.100;0.101;0.101;0.101;0.101;0.101;0.101;0.101;0.101;0.101;0.102; ...
0.102;0.102;0.102;0.102;0.102;0.102;0.102;0.103;0.103;0.103;0.103; ...
0.103;0.103;0.103;0.104;0.104;0.105;0.105;0.105;0.105;0.105;0.106; ...
0.106;0.106;0.106;0.107;0.107;0.107;0.108;0.108;0.108;0.108;0.108; ...
0.108;0.108;0.109;0.109;0.109;0.109;0.109;0.109;0.109;0.110;0.110; ...
0.110;0.110;0.110;0.110;0.110;0.110;0.110;0.110;0.110;0.111;0.111; ...
0.111;0.111;0.112;0.112;0.112;0.112;0.112;0.114;0.114;0.114;0.115; ...
0.115;0.115;0.117;0.119;0.119;0.119;0.119;0.120;0.120;0.120;0.121; ...
0.122;0.122;0.123;0.123;0.125;0.125;0.128; 0.129;0.129;0.129;0.130;0.131];
[H1, P1, CRITVAL1] = bootmode (x,1,2000);
% Repeat function call systematically increasing the number of modes (M) by
% 1, until the null hypothesis is accepted (i.e. H0 = 0)
[H2, P2, CRITVAL2] = bootmode (x,2,2000);
sprintf ('Summary of results:\n')
sprintf (cat (2, 'H1 is %u with p = %.3g so reject the null hypothesis', ...
'that there is 1 mode\n'), H1, P1)
sprintf (cat (2, 'H2 is %u with p = %.3g so accept the null hypothesis', ...
' that there are 2 modes\n'), H2, P2)
% Please be patient, these calculations take a while...
Produces the following output
ans = Summary of results: ans = H1 is 1 with p = 0.001 so reject the null hypothesisthat there is 1 mode ans = H2 is 0 with p = 0.35 so accept the null hypothesis that there are 2 modes
Package: statistics-resampling