function R = crr(s,varargin);
% CRR - Confidence Region Radius
%
%   CRR(S) computes the radius of a confidence interval, circle, or
%   sphere with 95% probability given S, a vector of standard deviations
%   from a multivariate normal distribution.
%
%   CRR(S,P) computes the confidence region radius with probability P.
%
%   CRR(S,P,M) performs a bootstrap validation with M normally
%   distributed random samples of size 1e6.
%
%   CRR(S,P,[M,N]) performs a bootstrap validation with M normally
%   distributed random samples of size N.
%
% EXAMPLES:
%
%   % 50% confidence interval
%   R=crr(1,0.5) % = 0.6745
%
%   % 95% confidence circle with bootstrap validation
%   R=crr([1,2],0.95,[10,1e5]) % = 4.0717
%
%   % 95% confidence circle from covariance matrix
%   C=[1 0.5; 0.5 2]; % = covariance matrix
%   R=crr(sqrt(eig(C))) % = 3.0915
%
%   % 90% confidence sphere
%   R=crr([1,1,1],0.9) % = 2.5003
%
%   % 95% confidence sphere from covariance matrix
%   C=[1 0.5 0.25; 0.5 2 -0.25; 0.25 -0.25 3]; % = covariance matrix
%   R=crr(sqrt(eig(C))) % = 4.0901
%
% REFERENCES:
%
%   Kleder, Michael (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, Ray (2004). "Joint normal distribution integral over the unit
%     sphere." The Math Forum.
%     (http://mathforum.org/kb/message.jspa?messageID=1552353&tstart=0)
%
% See also crr.pdf.

% Copyright(c)2006 Version 1.00
%   Tom Davis (tdavis@eng.usf.edu)
%   Michael Kleder (mkleder@hotmail.com)

msg='requires length(S) = 1, 2, or 3';
if nargin<1, error(msg), end
if length(s)>3 || length(s)<1, error(msg), end
smax=max(s);
if smax>=1e7*min(s), error('requires max(S)/min(S) < 1e7'), end
s=sort(s)/smax;
if nargin>1, p=varargin{1}; else, p=0.95; end
if isempty(p) || p<eps || p>1-eps, error('requires 0 < P < 1'), end
options=optimset('display','off','tolx',eps);
lowend=eps; highend=5;

n=length(s);
if n==1
  R=sqrt(2)*erfinv(p);
elseif n==2 & abs(diff(s))<eps
  R=sqrt(-2*log(1-p));
else
  warning('off','MATLAB:quadl:ImproperFcnValue')
  try
    R=fzero(@difference,[lowend,highend],options,s,p,1e-9);
  catch
    s=flipud(s(:));
    R=abs(fzero(@difference,lowend,options,s,p,1e-12));
  end
  warning('on','MATLAB:quadl:ImproperFcnValue')
end

if nargin>2 && ~isempty(varargin{2})
    bootstrap(R,s,varargin{2})
end

R=R*smax;

%--------------------------------------------------------------------
function d=difference(R,s,p,tol)
if length(s)==3 & abs(diff(s))<eps
  a=R/sqrt(2);
  cdf=erf(a)-2*a*exp(-a^2)/sqrt(pi);
elseif length(s)==3 && ...
  ((abs(s(1)-s(2))<eps && s(1)>s(3)) || ...
   (abs(s(2)-s(3))<eps && s(3)>s(1)))
  if s(3)>s(1), s=flipud(s); end
  a=R/sqrt(2);
  c=sqrt(s(1)^2-s(3)^2);
  cdf=erf(a/s(3))-s(1)/c*exp(-(a/s(1))^2)*erf(c*a/(s(1)*s(3)));
else
  cdf=quadl(@pdf,0,R,tol,0,R,s)/(s(1)*s(2));
end
d=p-cdf;

%--------------------------------------------------------------------
function f=pdf(r,R,s)
a= 0.25*(1/s(1)^2-1/s(2)^2);
b=-0.25*(1/s(1)^2+1/s(2)^2);
r2=r.^2;
f=r.*exp((b+abs(a))*r2);
if abs(a)>eps
  % Io(z,scale) = exp[-|Re(z)|] Io(z)
  f=f.*besseli(0,a*r2,1);
end
if length(s)==3
  f=f.*erf(sqrt(0.5*(R^2-r2))/s(3));
end

%--------------------------------------------------------------------
function bootstrap(R,s,mn)
disp('Running test of CRR.')
m=fix(mn(1));
if length(mn)>1, n=fix(mn(2)); else, n=1e6; end
j=length(s);
p=zeros(m,1);
for i=m:-1:1
  fprintf('%g ',i);
  x=randn(n,1)*s(1);
  if j>1, y=randn(n,1)*s(2); else, y=zeros(n,1); end
  if j>2, z=randn(n,1)*s(3); else, z=zeros(n,1); end
  r=sqrt(sum([x y z].^2,2));
  p(i)=sum(r<R)/length(r);
end
disp(' ');
disp([num2str(mean(p)*100),...
  '% of random points were located within the computed radius.'])