***--------------------------------------------------------------------***;
** Robustify data and Mardia's Skewness and Kurtosis;
**Technical details are provided in 
Yuan, Chan & Bentler (BJMSP 2000, 53, 31–50), 
and also in Yuan & Gomer (BJMSP 2020);
***--------------------------------------------------------------------***;
data Mardia;
*n=88;
infile 'd:\papers\2020\robustMethod\Mardia.dat';
input v1-v5; *5 variables for Mardia's data;
run;

options ls=80 nodate nonumber;
title;
proc iml;
reset noname;
*-----------------------------;
start order(n, m, t_stat);
*the order is from large to small according to the elements of the m column;
  do i=1 to n-1;
        do j=i+1 to n;
                if t_stat[i,m]<t_stat[j,m] then do;
                        tt=t_stat[i,];
                        t_stat[i,]=t_stat[j,];
                        t_stat[j,]=tt;
                end;
        end;
  end;
finish;

*-----------------------------;
*computing the sample skewness and 
kurtosis with unidimensional data;
*-----------------------------;
start sk(x,skew,kurt);
n=nrow(x);
x_bar=sum(x)/n;
m_x2=0;
m_x3=0;
m_x4=0;
x_c=x-j(n,1,x_bar);
x_c2=x_c#x_c;
m_x2=sum(x_c2)/n;
m_x3=sum(x_c2#x_c)/n;
m_x4=sum(x_c2#x_c2)/n;
skew=m_x3/(m_x2*sqrt(m_x2));
kurt=m_x4/(m_x2*m_x2);
finish;

*----------------------------------------------------;
*Computing Mardia's multivariate skewness and kurtosis;
*----------------------------------------------------;

start Mardia(x,mskew0,mkurt0,mkurt_p,df);
N=nrow(x);
p=ncol(x);

x_t=x`;
x_bar=x_t*j(N,1,1)/N;
Scov=(x_t*x-N*x_bar*x_bar`)/N;
mkurt=0.0;
scov_in=inv(Scov);

do i=1 to N;
        x_i=x[i,]`;
        x_i0=x_i-x_bar;
		x_it0=x_i0`;
        d_i2=x_it0*scov_in*x_i0;
        mkurt=mkurt+d_i2*d_i2;
end;
mkurt=mkurt/N;
mkurt_p=mkurt/(p*(p+2)); *relative Mardia's kurtosis;

mkurt0=(mkurt-p*(p+2))/sqrt(8*p*(p+2)/N); *standardized Mardia's kurtosis;

mskew=0.0;
do i=1 to N;
        x_i=x[i,]`;
        x_i0=x_i-x_bar;
		x_it0=x_i0`;
        do j=1 to N;
                x_j=x[j,]`;
                x_j0=x_j-x_bar;

                d_ij2=x_it0*scov_in*x_j0;
                mskew=mskew+d_ij2*d_ij2*d_ij2;
        end;
end;
mskew=mskew/(N*N);

mskew0=N*mskew/6; *standardized Mardia's skewness;
df=p*(p+1)*(p+2)/6; *nominal degrees of freedom for standardized Mardia's skewness;

finish;


***--------------------------------------------------------------------***;
*** Robust transformation with Huber's type of weight;
***--------------------------------------------------------------------***;
start robmusig(varphi,x,hmu_t,hsigma, xr, err);
ep=.00001;
N=nrow(x);
p=ncol(x);
prob=1-varphi;
chip2=cinv(prob,p);
kappa=sqrt(chip2);
tau=( p*probchi(chip2,p+2)+ chip2*(1-prob) )/p;
*print "kappa, tau=" kappa tau;

xr=j(N,p,0);

**initial values;
sum_xt=j(1,p,0); sum_xx=j(p,p,0);
do i=1 to N;
        x_it=x[i,];
		x_i=x_it`;
        sum_xt=sum_xt+x_it;
        sum_xx=sum_xx+x_i*x_it;
end;
mu_t0=sum_xt/N;
sigma_0=(sum_xx-N*(mu_t0`)*mu_t0)/(N-1);

delta=.01;
do jj=1 to 500 while (delta>ep);
	sig_in=inv(sigma_0);
	sum_w1=0.0;
	mu_t=j(1,p,0);
	sigma=j(p,p,0);
	do i=1 to N;
		x_it=x[i,];
  		x_it0=x_it-mu_t0;
  		x_i0=x_it0`;
  		sigxi_0=sig_in*x_i0;
  		d_i2=x_it0*sigxi_0;
  		d_i=sqrt(d_i2);
**Huber weight functions;
  		if d_i<=kappa then do;
     		w_i1=1.0;
  		end;
  		else do;
     		w_i1=kappa/d_i;
  		end;
  		w_i2=w_i1*w_i1;
  		sum_w1=sum_w1+w_i1;
  		mu_t=mu_t+w_i1*x_it;
  		sigma=sigma+w_i2*x_i0*x_it0;
	*i;
	end;
	mu_t1=mu_t/sum_w1;
	sigma_1=sigma/(N*tau);
	delta=max(abs(mu_t1-mu_t0), abs(sigma_1-sigma_0));
	mu_t0=mu_t1;
	sigma_0=sigma_1;
end;
if jj<499 then do;
	err=0;
	hmu_t=mu_t1;
    hsigma=sigma_1;
	*print "robust mean=" mu1;
	*print "robust covariance=" sigma1;
	dw_vec=j(n,4,0);
	do i=1 to N;
		x_it=x[i,];
  		x_it0=x_it-mu_t0;
  		x_i0=x_it0`;
  		sigxi_0=sig_in*x_i0;
  		d_i2=x_it0*sigxi_0;
  		d_i=sqrt(d_i2);
**Huber weight functions;
  		if d_i<=kappa then do;
     		w_i1=1.0;
  		end;
  		else do;
     		w_i1=kappa/d_i;
  		end;
  		w_i2=w_i1*w_i1/tau;
		xr[i,]=sqrt(w_i2)*x_it0;
		dw_vec[i,]=i||d_i||w_i1||w_i2;
	end;
	mean_xr=j(1,N,1)*xr/N;
	xr=xr+j(N,1,1)*(hmu_t-mean_xr);
	run order(n, 2, dw_vec); *order the case according to M-distance from large to small;
	id=(1:n)`;
	print "order, i, d_i, w_i1, w_i2";
	print id dw_vec [format=f6.3];
end;
else do;
	err=1;
	print "The IRLS method of Huber-M does not converge within 500 iterations";
end;

finish;


***--------------------------------------------------------------------***;
*** The following is the main program;
***--------------------------------------------------------------------***;
start main;

print "--------------------------------------------------------------------";
print "Mardia's open-closed book data";
print "---------------------------------------------------------";
use Mardia;
read all var _num_ into xmat;
* Mechanics, Vectors, Algebra, Analysis, Statistics;
*V1, V2 are  closed book, V3 V4 V5 are open book;
close Mardia;
x=xmat;
n=nrow(xmat); 
p=ncol(xmat); print "n, p=" n p;

print "---------------------------------------";
*marginal skewness and kurtosis;
skew=j(1,p,0);
kurt=j(1,p,0);
skew_r=j(1,p,0);
kurt_r=j(1,p,0);
do j=1 to p;
	x_j=x[,j];
	run sk(x_j,skew_j,kurt_j);
	skew[j]=skew_j;
	kurt[j]=kurt_j;
end;
print " Mechanics, Vectors, Algebra, Analysis, Statistics";
print "univariate skewness and kurtosis of raw data";
v_id=1:p;
print v_id;
print "skewness=" skew[format=f6.3];
print "kurtosis=" kurt[format=f6.3];
print "---------------------------------------";

run Mardia(x,mskew_0,mkurt_0,mkurt_p,df);
print "standardized Mardia's measures of skewness and kurtosis with raw data";
	print "df for mskew_0=" df;
	print "mskew_0=" mskew_0 [format=f10.3];
	print "mkurt_0=" mkurt_0 [format=f8.3];
	print "relative multivariate kurtosis=" mkurt_p [format=f8.3];

print "-----------------------------------------------------------";

varphi=0.05; *Tuning parameter varphi subject to change;
print "Tuning parameter varphi with Huber-weights =" varphi;
run robmusig(varphi,x,hmu_t,hsigma, xr, err);
if err=0 then do;
	print "robust means \hat\mu";
	print hmu_t;
	print "robust covariance matrix hat\Sigma=";
	print hsigma;
	do j=1 to p;
		xr_j=xr[,j];
		run sk(xr_j,skew_rj,kurt_rj);
		skew_r[j]=skew_rj;
		kurt_r[j]=kurt_rj;
	end;
print "univariate skewness and kurtosis of robustified data";
print v_id;
print "skewness=" skew_r[format=f6.3];
print "kurtosis=" kurt_r[format=f6.3];
	
print "---------------------------------------";

	run Mardia(xr,mskew_r,mkurt_r,mkurt_rp,df);
	print "standardized Mardia's measures of skewness and kurtosis with robustified data";

	print "df for mskew_r=" df;
	print "mskew_r=" mskew_r [format=f10.3];
	print "mkurt_r=" mkurt_r [format=f8.3];
	print "relative mkurt=" mkurt_rp [format=f8.3];
	print "-----------------------------------------------------------";
	print "Robustified data";
	print xr;
	print "-----------------------------------------------------------";
end;


finish;
run main;