Help with debugging an error: Attempted to access R(1); index out of bounds because numel(R)=0???

• Hi, I'm trying to run a bootstrap code that draws a random subsample from a dataset to use for bootstrapping. But now I get this error sequence:*

Attempted to access R(1); index out of bounds because numel(R)=0.

Error in regress (line 76)
p = sum(abs(diag(R)) > max(n,ncolX)*eps(R(1)));

Error in @(bootr)regress(Yfit+bootr,[ones(size(x,1),1),x])

Error in bootstrp (line 165)
    bootstat = feval(bootfun,bootargs{:});

Error in RandomSubsampleAnalysis (line 46)
[bootstat2]=bootstrp(rep,@(bootr)regress(Yfit+bootr,[ones(size(x,1),1) x]),resid);

Caused by:
    Unable to evaluate BOOTFUN with the supplied arguments.

The funny thing is that it worked fine a couple of hours ago.

my code looks like this:

clc;
warning off all;
replications=1; % number of simulations %

rep=100; % bootstrap replications %

rt1_s=[];
rt2_s=[];
OLS_b_s=[];
OLS_std_b_s=[];
Boot1_b_s=[];
Boot2_b_s=[];
Boot1_std_b_s=[];
Boot2_std_b_s=[];

for i=1:replications

% set subsample size %
ss=100; % cannot exceed full sample size %

load('data.mat');
n=length(data);
index=randsample(n,ss);
data1=data(index,:);
x=data1(:,2:end);
y=data1(:,1);

stats=regstats(y,x,'linear');
t=stats.tstat;

% case resampling

Y=[y ones(size(x,1),1) x];

[bootstat1]=bootstrp(rep,@(x)regress(x(:,1),x(:,2:end)),Y);

subplot(2,2,1), hist(bootstat1(:,1),100)
subplot(2,2,2), hist(bootstat1(:,2),100)

% model-based resampling

Yfit=[ones(size(x,1),1) x]*t.beta;
resid=y-Yfit;
[bootstat2]=bootstrp(rep,@(bootr)regress(Yfit+bootr,[ones(size(x,1),1) x]),resid);

subplot(2,2,3), hist(bootstat2(:,1),100)
subplot(2,2,4), hist(bootstat2(:,2),100)

beta = t.beta; % coefficients - OLS
se = t.se; % standard errors - OLS

Boot1_b=mean(bootstat1)'; % coefficients - case resampling
Boot1_std_b=sqrt(var(bootstat1))'; % standard errors - case resampling
rt1=abs((Boot1_b-t.beta)./Boot1_std_b);

Boot2_b=mean(bootstat2)'; % coefficients - model-based resampling
Boot2_std_b=sqrt(var(bootstat2))'; % standard errors - model-based resampling
rt2=abs((Boot2_b-t.beta)./Boot2_std_b);

rt1_s=[rt1_s;rt1'];
rt2_s=[rt2_s;rt2'];
OLS_b_s=[OLS_b_s;t.beta'];
OLS_std_b_s=[OLS_std_b_s;t.se'];
Boot1_b_s=[Boot1_b_s;Boot1_b'];
Boot2_b_s=[Boot2_b_s;Boot2_b'];
Boot1_std_b_s=[Boot1_std_b_s;Boot1_std_b'];
Boot2_std_b_s=[Boot2_std_b_s;Boot2_std_b'];

end

OLS_b=mean(OLS_b_s)'
OLS_std_b=mean(OLS_std_b_s)'

Boot1_b=mean(Boot1_b_s)'
Boot2_b=mean(Boot2_b_s)'

Boot1_std_b=mean(Boot1_std_b_s)'
Boot2_std_b=mean(Boot2_std_b_s)'

rt1=mean(rt1_s)'
rt2=mean(rt2_s)'

can someone help please?

Answer 1

[Stephen]

Are bootstrp and bootr functions you've defined? It would be really difficult to debug without seeing those. If the code used to work, and you haven't changed the functions, then I would suggest you should check the "Last Modified" date on the file 'data.mat' to see if that has changed since the last time it worked.

Good Luck,

Stephen

[IAn]

bootstrp is actually in the functions directory of 'bootstrap statistics'.

function [bootstat, bootsam] = bootstrp(nboot,bootfun,varargin)

%BOOTSTRP Bootstrap statistics.

% BOOTSTAT = BOOTSTRP(NBOOT,BOOTFUN,D1,...) draws NBOOT bootstrap data

% samples, computes statistics on each sample using the function BOOTFUN,

% and returns the results in the matrix BOOTSTAT. NBOOT must be a

% positive integer. BOOTFUN is a function handle specified with @.

% Each row of BOOTSTAT contains the results of applying BOOTFUN to one

% bootstrap sample. If BOOTFUN returns a matrix or array, then this

% output is converted to a row vector for storage in BOOTSTAT.

%

% The third and later input arguments (D1,...) are data (scalars,

% column vectors, or matrices) that are used to create inputs to BOOTFUN.

% BOOTSTRP creates each bootstrap sample by sampling with replacement

% from the rows of the non-scalar data arguments (these must have the

% same number of rows). Scalar data are passed to BOOTFUN unchanged.

%

% [BOOTSTAT,BOOTSAM] = BOOTSTRP(...) returns BOOTSAM, an N-by-B matrix of

% indices into the rows of the extra arguments, where N is the number of

% rows in non-scalar input arguments to BOOTSTRP and B is the number of

% generated bootstrap replicates. To get the output samples BOOTSAM

% without applying a function, set BOOTFUN to empty ([]).

%

% BOOTSTAT = BOOTSTRP(..., 'PARAM1',val1, 'PARAM2',val2, ...) specifies

% optional parameter name/value pairs to control how BOOTSTRP performs

% computations. Parameter names/values may only appear after the data

% arguments used as inputs to BOOTFUN. Parameters are:

%

% 'Weights' Observation weights. This must be a vector of

% non-negative numbers with at least one positive

% element. The number of elements in WEIGHTS must be

% equal to the number of rows in non-scalar input

% arguments to BOOTSTRP. To obtain one bootstrap

% replicate, BOOTSTRP samples N out of N with

% replacement using these weights as multinomial

% sampling probabilities.

%

% 'Options' A structure that contains options specifying whether to

% compute bootstrap iterations in parallel, and specifying

% how to use random numbers during the bootstrap sampling.

% This argument can be created by a call to STATSET.

% BOOTSTRP uses the following fields:

% 'UseParallel'

% 'UseSubstreams'

% 'Streams'

% For information on these fields see PARALLELSTATS.

% NOTE: if 'UseParallel' is 'always' and 'UseSubstreams'

% is 'never', then the length of Streams must equal the number

% of processors used by BOOTSTRP. There are two possibilities.

% If a MATLAB pool is open, then Streams is the same length as

% the size of the MATLAB pool. If a MATLAB pool is not open,

% then Streams must supply a single random number stream.

%

% The returned values from BOOTFUN are typically numeric values, but can

% be any values capable of being stored in an array. If the function

% returns objects (such as probability distributions) that cannot be

% stored in ordinary arrays, it can enclose those objects in a cell

% array.

%

% Examples:

%

% Compute a sample of 100 bootstrapped means of random samples taken from

% the vector Y, and plot an estimate of the density of these bootstrapped

% means:

% y = exprnd(5,100,1);

% m = bootstrp(100, @mean, y);

% [fi,xi] = ksdensity(m);

% plot(xi,fi);

%

% Compute a sample of 100 bootstrapped means and standard deviations of

% random samples taken from the vector Y, and plot the bootstrap estimate

% pairs:

% y = exprnd(5,100,1);

% stats = bootstrp(100, @(x) [mean(x) std(x)], y);

% plot(stats(:,1),stats(:,2),'o')

%

% Estimate the standard errors for a coefficient vector in a linear

% regression by bootstrapping residuals:

% load hald ingredients heat

% x = [ones(size(heat)), ingredients];

% y = heat;

% b = regress(y,x);

% yfit = x*b;

% resid = y - yfit;

% se = std(bootstrp(1000, @(bootr) regress(yfit+bootr,x), resid));

%

% Bootstrap a correlation coefficient standard error:

% load lawdata gpa lsat

% se = std(bootstrp(1000,@corr,gpa,lsat));

%

% Compute a sample of 100 bootstrapped means and standard deviations of

% random samples taken from the vector Y. Compute the bootstrap iterations

% in parallel (this only works with the Parallel Computing Toolbox).

% Plot the bootstrap estimate pairs:

% y = exprnd(5,100,1);

% matlabpool open; % only works with the Parallel Computing Toolbox

% opt = statset('UseParallel','always');

% stats = bootstrp(100, @(x) [mean(x) std(x)], y, 'Options', opt);

% plot(stats(:,1),stats(:,2),'o')

%

% Use bootstrapping to fit normal distributions to each sample, and

% collect the fitted distributions. Later examine values such as the

% interquartile range computed from those distributions:

%

% load census

% probs = bootstrp(100, @(x) {fitdist(x,'normal')}, pop);

% iqrs = cellfun(@iqr,probs);

% hist(iqrs)

%

% See also RANDOM, RANDSTREAM, STATSET, STATGET, PARALLELSTATS,

% HIST, KSDENSITY.

% Reference:

% Efron, Bradley, & Tibshirani, Robert, J.

% "An Introduction to the Bootstrap",

% Chapman and Hall, New York. 1993.

% Copyright 1993-2011 The MathWorks, Inc.

% $Revision: 1.1.8.4 $ $Date: 2011/05/09 01:24:12 $

% Sanity check the two initial arguments

if nargin<2

error(message('stats:bootstrp:TooFewInputs'));

end

if nboot<=0 || nboot~=round(nboot)

error(message('stats:bootstrp:BadNboot'))

end

% === Extract name-value pairs that are not arguments for bootfun ====

[weights, options, bootargs] = extractNameValuePairs(varargin{:});

% === Process the arguments to bootfun ===

[n,booteval] = bootEvalCommand(bootfun,bootargs{:});

% === Process the Options parameters ===

[useParallel, RNGscheme, poolsz] = ...

internal.stats.parallel.processParallelAndStreamOptions(options,true);

usePool = useParallel && poolsz>0;

% === Set up the bootstrp sampling function ===

[myrand,randargs] = defineRNGcall(RNGscheme, usePool, n, weights);

% === Begin actual processing ===

% Preallocate index matrix of bootstrap samples if requested.

haveallsamples = (nargout>=2);

% If no bootfun was supplied, we still generate a matrix of bootstrap samples,

% if two output arguments were supplied. This behavior is specified in

% the help text. We also return a dimensioned but empty matrix of zeros

% for the bootstrap statistics.

%

if isempty(bootfun)

bootstat = zeros(nboot,0);

if haveallsamples

bootsam = internal.stats.parallel.smartForSliceout( ...

nboot, @loopBodyEmptyBootfun, useParallel, RNGscheme);

end

return

end

% Sanity check bootfun call and determine dimension and type of result

try

% Get result of bootfun on actual data, force to a row.

bootstat = feval(bootfun,bootargs{:});

bootstat = bootstat(:)';

catch ME

m = message('stats:bootstrp:BadBootFun');

MEboot = MException(m.Identifier,'%s',getString(m));

ME = addCause(ME,MEboot);

rethrow(ME);

end

% Initialize an array to contain the results of all the bootstrap

% calculations, preserving the output type

bootstat(nboot,1:numel(bootstat)) = bootstat;

% === Do bootfun - nboot times. ===

if haveallsamples

[bootstat,bootsam] = internal.stats.parallel.smartForSliceout( ...

nboot, @loopBody, useParallel, RNGscheme);

else

bootstat = internal.stats.parallel.smartForSliceout( ...

nboot, @loopBody, useParallel, RNGscheme);

end

% ---- Nested functions ----

function onesample = loopBodyEmptyBootfun(~,~)

onesample = myrand(randargs{:});

end

function [onebootstat,onebootsam] = loopBody(~,~)

onesample = myrand(randargs{:});

tmp = booteval(onesample);

onebootstat = (tmp(:))';

if nargout > 1

onebootsam = onesample;

end

end

end % of bootstrp

function [myrand,randargs] = defineRNGcall(RNGscheme,usePool,n,weights)

uuid = RNGscheme.uuid;

streams = RNGscheme.streams;

useSubstreams = RNGscheme.useSubstreams;

if isempty(streams)

if isempty(weights)

myrand = @randi;

randargs = {n,n,1};

else

myrand = @randsample;

randargs = {n,n,true,weights};

end

else

% a stream or streams were supplied in the command line

S = streams{1};

if isempty(weights)

myrand = @randi;

if ~usePool || useSubstreams

% a single stream is in use throughout

randargs = {S,n,n,1};

else

%- We defer the stream assignment to within the loop iteration,

% the stream to be used by the worker is not known now

myrand = @(streams,useSubstreams,usePool) ...

randi(internal.stats.parallel.workerGetValue(uuid),n,n,1);

randargs = {streams,useSubstreams,usePool};

end

else

if ~usePool || useSubstreams

% a single stream is in use throughout

myrand = @randsample;

randargs = {S,n,n,true,weights};

else

%- We defer the stream assignment to within the loop iteration,

% the stream to be used by the worker is not known now

myrand = @(streams,useSubstreams,usePool) ...

randsample(internal.stats.parallel.workerGetValue(uuid),n,n,true,weights);

randargs = {streams,useSubstreams,usePool};

end

end

end

end %-defineRNGcall

function tmp = generalEval(onesample,bootfun,la,scalard,varargin)

db = cell(la,1);

for k = 1:la

if scalard(k) == 0

db{k} = varargin{k}(onesample,:);

else

db{k} = varargin{k};

end

end

tmp = feval(bootfun,db{:});

end %-generalEval

function [n,anonfun] = bootEvalCommand(bootfun,varargin)

% === Process the arguments to bootfun ===

% Initialize matrix to identify scalar arguments to bootfun.

la = length(varargin);

if la == 0

error(message('stats:bootstrp:NotEnoughBootfunArgs'));

end

scalard = zeros(la,1);

% find out the size information in varargin.

n = 1;

for k = 1:la

[row,col] = size(varargin{k});

if max(row,col) == 1

scalard(k) = 1;

end

if row == 1 && col ~= 1

row = col;

varargin{k} = varargin{k}(:);

end

if n>1 && row>1 && row~=n

error(message('stats:bootstrp:SizeMismatchBootfunArgs'));

end

n = max(n,row);

end

if n<2

error(message('stats:bootstrp:OnlyScalarBootfunArgs'));

end

% === Define anonymous function to evaluate bootfun with the supplied arguments ===

if la==1 && ~any(scalard)

X1 = varargin{1};

anonfun = @(onesample) feval(bootfun,X1(onesample,:));

elseif la==2 && ~any(scalard)

X1 = varargin{1};

X2 = varargin{2};

anonfun = @(onesample) feval(bootfun,X1(onesample,:),X2(onesample,:));

else

anonfun = @(onesample) generalEval(onesample,bootfun,la,scalard,varargin{:});

end

end %-bootEval

function [weights,options,bootargs] = extractNameValuePairs(varargin)

% Scan the argument list until we run into a string. Everything to the

% right of the string is considered arguments not for bootfun.

weights = [];

options = statset('bootstrp');

defSpecialArgs = {'weights' 'options'};

defSpecialValues = {weights options};

if length(varargin)>1

screen = @(arg) (ischar(arg) && size(arg,1)==1);

isspecial = cellfun(screen,varargin);

newOptions = [];

if any(isspecial)

firstspecial = find(isspecial,1,'first');

specialArgs = varargin(firstspecial:end);

varargin = varargin(1:firstspecial-1);

[weights,newOptions] ...

= internal.stats.parseArgs(defSpecialArgs,defSpecialValues,specialArgs{:});

end

% Check parallel options

if ~isempty(newOptions)

if ~isstruct(newOptions)

error(message('stats:bootstrp:BadOption'))

end

try

options = statset(options, newOptions);

catch ME

m = message('stats:bootstrp:BadOption');

newME = MException(m.Identifier,'%s',getString(m));

newME = addCause(newME,ME);

throw(newME);

end

end

end

bootargs = varargin;

end %-extractNameValuePairs

I doubt the problem is with my data.mat file because it is simply a matrix of values for my regressors? it almost seems like the error is saying that the problem is with the supplied bootstrp function that comes with matlab to begin with?

[Stephen]

Unfortunately, I don't have the statistics toolbox so I don't have thee bootr function either. Have you tried adding breakpoints to the script and stepping in to the function to find where the issue occurs and what the value of R is when the error occurs?

Stephen