function R = crr(s,varargin)
% CRR - Confidence Region Radius (R13)
%
% R = CRR(S) computes the radius of the mean-centered interval, circle,
% or sphere with 95% probability given S, which is either a vector of
% standard deviations or a covariance matrix from a multivariate normal
% distribution. If S is a real, symmetric, positive semidefinite
% matrix, CRR(S) is equivalent to CRR(SQRT(EIG(S))). Scalar S is treated
% as a standard deviation.
%
% R = CRR(S,P) computes the confidence region radius with probability P
% instead of the default, which is 0.95.
%
% R = CRR(S,P,TOL) uses a quadrature tolerance of TOL instead of the
% default, which is 1e-15. Larger values of TOL may result in fewer
% function evaluations and faster computation, but less accurate results.
% Use [] as a placeholder to obtain the default value of P.
%
% R = CRR(S,P,TOL,M) performs a bootstrap validation with M normally
% distributed random samples of size 1e6. Use [] as a placeholder to
% obtain the default value of TOL.
%
% R = CRR(S,P,TOL,[M N]) performs a bootstrap validation with M normally
% distributed random samples of size N.
%
% EXAMPLES:
%
% % Probable Error (PE)
% r=crr(1,0.5) % = 0.6745
%
% % 95% confidence interval
% r=crr(1) % = 1.9600
%
% % 90% confidence interval from data vector X
% r=crr(std(X),0.9)
%
% % Circular Error Probable (CEP)
% r=crr([1 2],0.5) % = 1.7408
%
% % 95% confidence circle with bootstrap validation
% r=crr([1 2],[],[],[10 1e5]) % = 4.0717
%
% % 90% confidence circle from covariance matrix V
% V=[1 1; 1 4];
% r=crr(V,0.9) % = 3.5263
%
% % 95% confidence circle from data matrix [X Y]
% r=crr(cov([X Y]))
%
% % Spherical Error Probable (SEP)
% r=crr([1 2 3],0.5) % = 3.1068
%
% % 95% confidence sphere from covariance matrix V
% V=[1 1 2; 1 4 -1; 2 -1 9];
% r=crr(V) % = 6.5817
%
% % 90% confidence sphere from data matrix [X Y Z]
% r=crr(cov([X Y Z]),0.9)
%
% % Inverse chi-square
% % chi2inv(p,n) = crr(eye(n),p)^2 = crr(ones(1,n),p)^2; n=1,2,3
% r2=crr(1)^2 % = 3.8415 = chi2inv(0.95,1)
% r2=crr([1 1])^2 % = 5.9915 = chi2inv(0.95,2)
% r2=crr([1 1 1])^2 % = 7.8147 = chi2inv(0.95,3)
%
% REFERENCES:
%
% Abramowitz, M. and Stegun, I. A., eds. (A&S 1972). Handbook of
% mathematical functions with formulas, graphs, and mathematical
% tables. National Bureau of Standards. 6.5.1, 7.1.23, 7.1.29, 9.6.12.
% (http://www.math.sfu.ca/~cbm/aands/frameindex.htm)
% Acklam, P. J. (2002). "Variance matrix checking: varchk."
% (http://home.online.no/~pjacklam/matlab/software/util/statutil/
% varchk.m)
% Altham, P. M. E. (2006). "Applied multivariate analysis, notes."
% (http://www.statslab.cam.ac.uk/~pat/AppMultNotes.pdf)
% 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)
% Godfrey, P. and Acklam, P. J. (2003). "Error function for complex
% inputs: erfz." Mathworks Central File Exchange.
% (http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?
% objectId=3574&objectType=file)
% Kleder, M. (2004). "An algorithm for converting covariance to
% spherical error probable." Mathworks Central File Exchange.
% (http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?
% objectId=5688&objectType=file)
% Koopman, R. (2004). "Joint normal distribution integral over the unit
% sphere." The Math Forum.
% (http://mathforum.org/kb/message.jspa?messageID=1552353&tstart=0)
% National Institute of Standards and Technology (NIST 2006).
% "The multivariate normal distribution."
% (http://www.itl.nist.gov/div898/handbook/pmc/section5/pmc542.htm)
% Sloane, N. J. A. (2006). "The on-line encyclopedia of integer
% sequences." AT&T Labs Research. A000984
% (http://www.research.att.com/~njas/sequences/A000984)
%
% See also crr.pdf.
% Copyright(c)2006-2008
% Tom Davis (tdavis@metzgerwillard.com)
% Michael Kleder (mkleder@hotmail.com)
%
% Last revision: 03/16/2008
if nargin<1 || isempty(s)
error('CRR 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 length(s)>3;
error('Covariance matrix cannot have more than three rows.')
elseif ~isreal(s)
error('Covariance matrix must be real.')
elseif any(any(s~=s.'))
error('Covariance matrix must be symmetric.')
end
s=eig(s);
if any(s<-eps)
error('Covariance matrix must be positive semidefinite.')
end
s(abs(s)<=eps)=0; s=sqrt(s);
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 length(s)>3;
error('Deviation vector cannot have more than three elements.')
elseif ~isreal(s)
error('Deviation vector must be real.')
elseif any(s<-eps)
error('Deviation vector must be nonnegative.')
end
s(abs(s)<=eps)=0;
end
if all(s==0), R=0; return, end
s=s(s~=0); smax=max(s);
if smax>=1e7*min(s)
error('Ratio of nonzero deviations must be less than 1e7.')
end
s=sort(s)/smax;
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('CRR requires 0 < P < 1.'), end
if nargin>2 && ~isempty(varargin{2}) && isa(varargin{2},'double') ...
&& numel(varargin{2})==1 && real(varargin{2})~=0
tol=real(varargin{2}); else tol=1e-15; % quadrature tolerance
end
if tol<0, tol=-tol; qflag=0; else qflag=1; end % quadrature flag
options=optimset('display','off','tolx',eps); % fzero tolerance
n=length(s);
if n==1 % crr.pdf (3)
R=sqrt(2)*erfinv(p);
elseif n==2 && abs(diff(s))<eps % crr.pdf (10)
R=sqrt(-2*log(1-p));
elseif n==3 && all(abs(diff(s))<eps) % crr.pdf (18)
R=sqrt(-2*log(1-p)); R0=1e20;
b=sqrt(pi)/2; count=0;
while abs(R-R0)>tol && count<50 % Newton's method
a=R/sqrt(2); R0=R; count=count+1;
R=R+(1+b*exp(a^2)*(p-erf(a))/a)/R;
end
if count==50, R=(R+R0)/2; end
else
R0=sqrt(2)*erfinv(p);
try % crr.pdf (7,9,15,16)
R=fzero(@difference,[R0 R0+1],options,s,p,tol,qflag);
catch % crr.pdf (7,15)
s=flipud(s(:));
R=abs(fzero(@difference,R0,options,s,p,tol,0));
end
end
if nargin>3 && ~isempty(varargin{3}) && isa(varargin{3},'double')
bootstrap(R,s,varargin{3})
end
R=R*smax;
%--------------------------------------------------------------------------
function f=difference(R,s,p,tol,qflag)
n=length(s);
if qflag && n==3 && any(abs(diff(s))<eps) % crr.pdf (16)
if abs(s(3)-s(2))<eps, s=[s(2) s(3) s(1)]; end
a=R/sqrt(2); b=a/s(3);
c=sqrt(1-(s(3)/s(1))^2);
F=erf(b)-exp(-(a/s(1))^2)*erfir(c*b)/c;
elseif qflag && n==2 && max(s)/min(s)<=10 && p<=0.995 % crr.pdf (9)
k=160; K=2*(0:k)+1;
a=0.25*(s(1)^-2-s(2)^-2);
b=0.25*(s(1)^-2+s(2)^-2);
% central binomial coefficients (cbc)
% C(2k,k) = (2k)!/(k!)^2 (Sloane A000984)
cbc=ones(1,k+1);
for n=1:k, cbc(n+1)=cbc(n)*(4-2/n); end
F=sum(cbc.*(0.5*a/b).^K.*gammainc(b*R^2,K))/(a*s(1)*s(2));
else % crr.pdf (7,15)
F=quadl(@integrand,0,R,tol,0,R,s)/(s(1)*s(2));
end
f=p-F;
%--------------------------------------------------------------------------
function f=integrand(r,R,s)
n=length(s);
a=abs(0.25*(s(1)^-2-s(2)^-2));
b=-0.25*(s(1)^-2+s(2)^-2);
r2=r.^2;
f=r.*exp((b+a)*r2); % crr.pdf (7,15)
if a>eps
% Io( z,1) = exp[-|Re(z)|] Io(z)
% Io(-z,1) = Io(z,1)
f=f.*besseli(0,a*r2,1); % crr.pdf (7,15)
end
if n==3
f=f.*erf(sqrt(0.5*(R^2-r2))/s(3)); % crr.pdf (15)
end
%--------------------------------------------------------------------------
function f=erfir(z)
% error function of pure imaginary or real argument;
% cf. erf(z) = gammainc(z^2,0.5)
y=imag(z);
if abs(y)<eps
f=erf(real(z));
elseif abs(y)<=8 % A&S 7.1.29
n=1:32;
f=i/pi*(sum(exp(-0.25*n.*n)./n.*sinh(n*y))*2+y);
elseif abs(y)<27 % A&S 7.1.23
m=193; s=1; y2=2*y^2;
for n=m:-2:1, s=1+n*s/y2; end
f=i*s*exp(y^2)/(y*sqrt(pi));
else
f=sign(y)*i*inf;
end
%--------------------------------------------------------------------------
function bootstrap(R,s,mn)
disp('Running test of CRR.')
mn=fix(abs(real(mn)));
m=max(mn(1),1);
if length(mn)>1, n=max(mn(2),1); else n=1e6; end
j=length(s);
p=zeros(m,1);
for i=m:-1:1
fprintf('%g ',i);
X=randn(n,j)*diag(s);
p(i)=sum(sum(X.^2,2)<R^2)/n;
end
disp(' ');
disp([num2str(mean(p)*100),'% of ',num2str(n*m),...
' random points were located within the computed radius.'])
