fastmcd.m

SVDD的工具箱

         line(ellip(:,1),ellip(:,2));
      end
      
   end
   
   if allsubset & choice < 6
      ask=0;
      choice=choice+1;
   elseif allsubset==1 & choice==6
      choice=7;
   elseif allsubset==2 & choice==6
      allsubset=0;
      ask=1;
      closeplot=1;
   end
   
 end
   
%----------------------------------------------------------------------------------------   

function beg(xlab,ylab,ord,x,y,nid,n,xmin,xmax,ymin,ymax)

% Creates a scatter plot.

scatter(x,y,3,'k')
      
xlabel(xlab);
ylabel(ylab);

xlim([xmin,xmax]);
ylim([ymin,ymax]);
box;
if nid   
   [ord,ind]=sort(ord);
   ind=ind(n-nid+1:n)';
	text(x(ind),y(ind),int2str(ind));
end

%----------------------------------------------------------------------------------------

function coord=ellipse(mean,covar)

% Determines the coordinates of some points that lie on the 97.5 % tolerance ellipse.

deter=covar(1,1)*covar(2,2)-covar(1,2)^2;
ylimit=sqrt(7.37776*covar(2,2));
y=-ylimit:0.005*ylimit:ylimit;
sqtdi=sqrt(deter*(ylimit^2-y.^2))/covar(2,2);
sqtdi([1,end])=0;
b=mean(1)+covar(1,2)/covar(2,2)*y;
x1=b-sqtdi;
x2=b+sqtdi;
y=mean(2)+y;
coord=[x1,x2([end:-1:1]);y,y([end:-1:1])]';

%-----------------------------------------------------------------------------------------

function quan=quanf(alpha,n,rk)

quan=floor(2*floor((n+rk+1)/2)-n+2*(n-floor((n+rk+1)/2))*alpha);

%-----------------------------------------------------------------------------------------

function rawconsfac=rawconsfactor(quan,n,p)

qalpha=qchisq(quan/n,p);
calphainvers=pgamma(qalpha/2,p/2+1)/(quan/n);
calpha=1/calphainvers;
rawconsfac=calpha;

%-----------------------------------------------------------------------------------------

function rewconsfac=rewconsfactor(weights,n,p)

if sum(weights)==n
   cdelta.rew=1;
else
   qdelta.rew=qchisq(sum(weights)/n,p);
   cdeltainvers.rew=pgamma(qdelta.rew/2,p/2+1)/(sum(weights)/n);
   cdelta.rew=1/cdeltainvers.rew;
end
rewconsfac=cdelta.rew;

%-----------------------------------------------------------------------------------------

function rawcorfac=rawcorfactor(p,n,alpha)

if p > 2 
   coeffqpkwad875=[-0.455179464070565,1.11192541278794,2;-0.294241208320834,1.09649329149811,3]';
   coeffqpkwad500=[-1.42764571687802,1.26263336932151,2;-1.06141115981725,1.28907991440387,3]';
   y1_500=1+(coeffqpkwad500(1,1)*1)/p^coeffqpkwad500(2,1);
   y2_500=1+(coeffqpkwad500(1,2)*1)/p^coeffqpkwad500(2,2);
   y1_875=1+(coeffqpkwad875(1,1)*1)/p^coeffqpkwad875(2,1);
   y2_875=1+(coeffqpkwad875(1,2)*1)/p^coeffqpkwad875(2,2);
   y1_500=log(1-y1_500);
	y2_500=log(1-y2_500);
   y_500=[y1_500;y2_500];
   A_500=[1,log(1/(coeffqpkwad500(3,1)*p^2));1,log(1/(coeffqpkwad500(3,2)*p^2))];
   coeffic_500=A_500\y_500;
   y1_875=log(1-y1_875);
	y2_875=log(1-y2_875);
   y_875=[y1_875;y2_875];
   A_875=[1,log(1/(coeffqpkwad875(3,1)*p^2));1,log(1/(coeffqpkwad875(3,2)*p^2))];
   coeffic_875=A_875\y_875;
   fp_500_n=1-(exp(coeffic_500(1))*1)/n^coeffic_500(2);
   fp_875_n=1-(exp(coeffic_875(1))*1)/n^coeffic_875(2);
else 
   if p == 2
      fp_500_n=1-(exp(0.673292623522027)*1)/n^0.691365864961895;
      fp_875_n=1-(exp(0.446537815635445)*1)/n^1.06690782995919;
   end   
   if p == 1 
      fp_500_n=1-(exp(0.262024211897096)*1)/n^0.604756680630497;
      fp_875_n=1-(exp(-0.351584646688712)*1)/n^1.01646567502486;
   end   
end   
if 0.5 <= alpha & alpha <= 0.875 
   fp_alpha_n=fp_500_n+(fp_875_n-fp_500_n)/0.375*(alpha-0.5);
end         
if 0.875 < alpha & alpha < 1
   fp_alpha_n=fp_875_n+(1-fp_875_n)/0.125*(alpha-0.875);
end         
rawcorfac=1/fp_alpha_n;

%-----------------------------------------------------------------------------------------

function rewcorfac=rewcorfactor(p,n,alpha)

if p > 2 
   coeffrewqpkwad875=[-0.544482443573914,1.25994483222292,2;-0.343791072183285,1.25159004257133,3]';
   coeffrewqpkwad500=[-1.02842572724793,1.67659883081926,2;-0.26800273450853,1.35968562893582,3]';
   y1_500=1+(coeffrewqpkwad500(1,1)*1)/p^coeffrewqpkwad500(2,1);
   y2_500=1+(coeffrewqpkwad500(1,2)*1)/p^coeffrewqpkwad500(2,2);
   y1_875=1+(coeffrewqpkwad875(1,1)*1)/p^coeffrewqpkwad875(2,1);
   y2_875=1+(coeffrewqpkwad875(1,2)*1)/p^coeffrewqpkwad875(2,2);
	y1_500=log(1-y1_500);
	y2_500=log(1-y2_500);
   y_500=[y1_500;y2_500];
   A_500=[1,log(1/(coeffrewqpkwad500(3,1)*p^2));1,log(1/(coeffrewqpkwad500(3,2)*p^2))];
   coeffic_500=A_500\y_500;
   y1_875=log(1-y1_875);
	y2_875=log(1-y2_875);
   y_875=[y1_875;y2_875];
   A_875=[1,log(1/(coeffrewqpkwad875(3,1)*p^2));1,log(1/(coeffrewqpkwad875(3,2)*p^2))];
   coeffic_875=A_875\y_875;
   fp_500_n=1-(exp(coeffic_500(1))*1)/n^coeffic_500(2);
   fp_875_n=1-(exp(coeffic_875(1))*1)/n^coeffic_875(2);
else 
   if p == 2
      fp_500_n=1-(exp(3.11101712909049)*1)/n^1.91401056721863;
      fp_875_n=1-(exp(0.79473550581058)*1)/n^1.10081930350091;
   end
   if p == 1 
      fp_500_n=1-(exp(1.11098143415027)*1)/n^1.5182890270453;
      fp_875_n=1-(exp(-0.66046776772861)*1)/n^0.88939595831888;
   end
end
if 0.5 <= alpha & alpha <= 0.875 
   fp_alpha_n=fp_500_n+(fp_875_n-fp_500_n)/0.375*(alpha-0.5);
end            
if 0.875 < alpha & alpha < 1 
   fp_alpha_n=fp_875_n+(1-fp_875_n)/0.125*(alpha-0.875);
end            
rewcorfac=1/fp_alpha_n;

%-----------------------------------------------------------------------------------------

function x = qchisq(p,a)
%QCHISQ   The chisquare inverse distribution function
%
%         x = qchisq(p,DegreesOfFreedom)

%        Anders Holtsberg, 18-11-93
%        Copyright (c) Anders Holtsberg

if any(any(abs(2*p-1)>1))
   error('A probability should be 0<=p<=1, please!')
end
if any(any(a<=0))
   error('DegreesOfFreedom is wrong')
end

x = qgamma(p,a*0.5)*2;

%-----------------------------------------------------------------------------------------

function x = qgamma(p,a)
%QGAMMA   The gamma inverse distribution function
%
%         x = qgamma(p,a)

%        Anders Holtsberg, 18-11-93
%        Copyright (c) Anders Holtsberg

if any(any(abs(2*p-1)>1))
   error('A probability should be 0<=p<=1, please!')
end
if any(any(a<=0))
   error('Parameter a is wrong')
end

x = max(a-1,0.1);
dx = 1;
while any(any(abs(dx)>256*eps*max(x,1)))
   dx = (pgamma(x,a) - p) ./ dgamma(x,a);
   x = x - dx;
   x = x + (dx - x) / 2 .* (x<0);
end

I0 = find(p==0);
x(I0) = zeros(size(I0));
I1 = find(p==1);
x(I1) = zeros(size(I0)) + Inf;

%-----------------------------------------------------------------------------------------

function f = dgamma(x,a)
%DGAMMA   The gamma density function
%
%         f = dgamma(x,a)

%       Anders Holtsberg, 18-11-93
%       Copyright (c) Anders Holtsberg

if any(any(a<=0))
   error('Parameter a is wrong')
end

f = x .^ (a-1) .* exp(-x) ./ gamma(a);
I0 = find(x<0);
f(I0) = zeros(size(I0));

%-----------------------------------------------------------------------------------------

function F = pgamma(x,a)
%PGAMMA   The gamma distribution function
%
%         F = pgamma(x,a)

%       Anders Holtsberg, 18-11-93
%       Copyright (c) Anders Holtsberg

if any(any(a<=0))
   error('Parameter a is wrong')
end

F = gammainc(x,a);
I0 = find(x<0);
F(I0) = zeros(size(I0));

%-----------------------------------------------------------------------------------------

function x = rchisq(n,a)
%RCHISQ   Random numbers from the chisquare distribution
%
%         x = rchisq(n,DegreesOfFreedom)

%        Anders Holtsberg, 18-11-93
%        Copyright (c) Anders Holtsberg

if any(any(a<=0))
   error('DegreesOfFreedom is wrong')
end

x = rgamma(n,a*0.5);


%-----------------------------------------------------------------------------------------

function x = rgamma(n,a)
%RGAMMA   Random numbers from the gamma distribution
%
%         x = rgamma(n,a)

%        Anders Holtsberg, 18-11-93
%        Copyright (c) Anders Holtsberg

if any(any(a<=0))
   error('Parameter a is wrong')
end

if size(n)==1
   n = [n 1];
end

x = qgamma(rand(n),a);

%-----------------------------------------------------------------------------------------

function normqqplot(x,y);
 
y = sort(y);

scatter(x,y,3,'k')

%-----------------------------------------------------------------------------------------

%function asvar=asvardiag(quan,n,p)
%
%alfa=quan/n;
%alfa=1-alfa;
%qalfa=qchisq(1-alfa,p);
%calfainvers=pgamma(qalfa/2,p/2+1);
%calfa=(1-alfa)/calfainvers;
%c2=-1/2*pgamma(qalfa/2,p/2+1);
%c3=-1/2*pgamma(qalfa/2,p/2+2);
%c4=3*c3;
%b1=(calfa*(c3-c4))/(1-alfa);
%b2=1/2+(calfa/(1-alfa))*(c3-((qalfa/p)*(c2+(1-alfa)/2)));
%asvar=(1-alfa)*b1^2*(alfa*((calfa*qalfa)/p-1)^2-1);
%asvar=asvar-2*c3*(calfa)^2*(3*(b1-p*b2)^2+(p+2)*b2*(2*b1-p*b2));
%asvar=asvar/(((1-alfa)*b1*(b1-p*b2))^2);
%
%----------------------------------------------------------------------------------------

function x = qf(p,a,b)
%QF       The F inverse distribution function
%
%         x = qf(p,df1,df2)

%        Anders Holtsberg, 18-11-93
%        Copyright (c) Anders Holtsberg

x = qbeta(p,a/2,b/2);
x = x.*b./((1-x).*a);

%----------------------------------------------------------------------------------------

function x = qbeta(p,a,b)
%QBETA    The beta inverse distribution function
%
%         x = qbeta(p,a,b)

%       Anders Holtsberg, 27-07-95
%       Copyright (c) Anders Holtsberg

if any(any((a<=0)|(b<=0)))
   error('Parameter a or b is nonpositive')
end
if any(any(abs(2*p-1)>1))
   error('A probability should be 0<=p<=1, please!')
end
b = min(b,100000);

x = a ./ (a+b);
dx = 1;
while any(any(abs(dx)>256*eps*max(x,1)))
   dx = (betainc(x,a,b) - p) ./ dbeta(x,a,b);
   x = x - dx;
   x = x + (dx - x) / 2 .* (x<0);
   x = x + (1 + (dx - x)) / 2 .* (x>1);
end

%-----------------------------------------------------------------------------------------

function d = dbeta(x,a,b)
%DBETA    The beta density function
%
%         f = dbeta(x,a,b)

%       Anders Holtsberg, 18-11-93
%       Copyright (c) Anders Holtsberg

if any(any((a<=0)|(b<=0)))
   error('Parameter a or b is nonpositive')
end

I = find((x<0)|(x>1));

d = x.^(a-1) .* (1-x).^(b-1) ./ beta(a,b);
d(I) = 0*I;

%----------------------------------------------------------------------------------------