function [r,e] = rcr(s,varargin)
% RCR - Rectangular Confidence Regions (R13)
%
% R = RCR(S) computes the semi-edge-length of the mean-centered hypercube
% with 95% probability given S, which is either a covariance matrix or a
% vector of standard deviations from a multivariate normal distribution.
% If S is a real, nonnegative vector, RCR(S) is equivalent to
% RCR(DIAG(S.^2)). Scalar S is treated as a standard deviation.
%
% R = RCR(S,P) computes the semi-edge-length of the hypercube with
% probability P instead of the default, which is 0.95. R is the two-
% tailed, equicoordinate quantile corresponding to P. The hypercube
% edge-length is 2*R.
%
% R = RCR(S,P,NP) uses NP quadrature points instead of the default, which
% is 2^11. Smaller values of NP result in faster computation, but may
% yield less accurate results. Use [] as a placeholder to obtain the
% default value of P.
%
% R = RCR(S,P,NP,M) performs a bootstrap validation with M normally
% distributed random samples of size 1e6. Use [] as a placeholder to
% obtain the default value of NP.
%
% R = RCR(S,P,NP,[M N]) performs a bootstrap validation with M normally
% distributed random samples of size N.
%
% [R,E] = RCR(S,...) returns an error estimate E.
%
% EXAMPLES:
%
% % 95% confidence interval
% r=rcr(1) % = 1.9600
%
% % 90% confidence interval from data vector X
% r=rcr(std(X),0.9)
%
% % 95% confidence square with bootstrap validation
% r=rcr([1 2],[],[],[10 1e5]) % = 3.9214
%
% % 90% confidence square from covariance matrix V
% V=[1 1; 1 4];
% r=rcr(V,0.9) % = 3.2939
%
% % 95% confidence square from data matrix [X Y]
% r=rcr(cov([X Y]))
%
% % 95% confidence cube from covariance matrix V
% V=[1 1 2; 1 4 -1; 2 -1 9];
% r=rcr(V) % = 5.9508
%
% % 90% confidence cube from data matrix [X Y Z]
% r=rcr(cov([X Y Z]),0.9)
%
% % 95% confidence hypercube with error estimate
% V=[4 1 2 1; 1 4 1 2; 2 1 2 1; 1 2 1 9];
% [r,e]=rcr(V,[],2^14) % = [6.0050, 1.3510e-005]
%
% REFERENCES:
%
% Acklam, P. J. (2002). "Variance matrix checking: varchk."
% (http://home.online.no/~pjacklam/matlab/software/util/statutil/
% varchk.m)
% Davis, T. (2006). "Rectangular confidence regions."
% Mathworks Central File Exchange.
% (http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?
% objectId=11627&objectType=file)
% Davis, T. and Kleder, M. (2006). "Confidence region radius."
% Mathworks Central File Exchange.
% (http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?
% objectId=10526&objectType=file)
% Genz, A. (2006). "Numerical computation of multivariate normal
% probabilities: qsimvnv."
% (http://www.math.wsu.edu/faculty/genz/software/matlab/qsimvnv.m)
% Genz, A. (1992). "Numerical computation of multivariate normal
% probabilities." J. Comp. Graph. Stat. 1, 141-149.
% (http://www.math.wsu.edu/faculty/genz/papers/mvn.pdf)
%
% See also qsimvnv.m and mvn.pdf.
% Copyright(c)2006-2008
% Tom Davis (tdavis@metzgerwillard.com)
%
% Last revision: 03/16/2008
if nargin<1 || isempty(s)
error('RCR requires at least one argument.')
end
if all(size(s)>1) % covariance matrix
if ~isa(s,'double')
error('Covariance matrix must be of class ''double''.')
elseif ndims(s)>2
error('Covariance matrix cannot have more than two dimensions.')
elseif size(s,1)~=size(s,2);
error('Covariance matrix must be square.')
elseif ~isreal(s)
error('Covariance matrix must be real.')
elseif any(any(s~=s.'))
error('Covariance matrix must be symmetric.')
end
ss=eig(s);
if any(ss<0)
error('Covariance matrix must be positive semi-definite.')
end
smax=max(sqrt(ss));
else % deviation vector
if ~isa(s,'double')
error('Deviation vector must be of class ''double''.')
elseif ndims(s)>2
error('Deviation vector cannot have more than two dimensions.')
elseif ~isreal(s)
error('Deviation vector must be real.')
elseif any(s<0)
error('Deviation vector must be nonnegative.')
else
smax=max(s); s=diag(s.^2);
end
end
if nargin>1 && ~isempty(varargin{1}) && isa(varargin{1},'double') ...
&& numel(varargin{1})==1
p=real(varargin{1}); else p=0.95; % probability
end
if p<eps || p>1-eps, error('RCR requires 0 < P < 1.'), end
if nargin>2 && ~isempty(varargin{2}) && isa(varargin{2},'double') ...
&& numel(varargin{2})==1 && real(varargin{2})~=0
np=fix(real(varargin{2})); else np=2^11; % number of points
end
options=optimset('display','off','tolx',eps); % fzero tolerance
R0=sqrt(2)*smax*erfinv(p);
r=abs(fzero(@difference,R0,options,s,p,np));
bound=abs(r)*ones(length(s),1);
[ignore,e]=qsimvnv(np,s,-bound,bound); %#ok ignore not used
if nargin>3 && ~isempty(varargin{3}) && isa(varargin{3},'double')
bootstrap(r,s,varargin{3})
end
%--------------------------------------------------------------------------
function f=difference(r,s,p,np)
bound=abs(r)*ones(length(s),1);
f=p-qsimvnv(np,s,-bound,bound);
%--------------------------------------------------------------------------
function bootstrap(r,s,mn)
disp('Running test of RCR.')
mn=fix(abs(real(mn)));
m=max(mn(1),1);
if length(mn)>1, n=max(mn(2),1); else n=1e6; end
[V,D]=eig(s);
j=length(s);
p=zeros(m,1);
for i=m:-1:1
fprintf('%g ',i);
X=randn(n,j)*sqrt(D)*V';
p(i)=sum(sum(abs(X)<r,2)==j)/n;
end
disp(' ');
disp([num2str(mean(p)*100),'% of ',num2str(n*m),...
' random points were located within the computed hypercube.'])
