MATLAB Answers

IAn
0

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

Asked by IAn
on 23 Mar 2012
  • 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?

  0 Comments

Sign in to comment.

Tags

Products

1 Answer

Answer by Stephen on 23 Mar 2012

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

  2 Comments

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?

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

Sign in to comment.