#** Technical details are provided in 
# Yuan, Chan & Bentler (BJMSP 2000, 53, 31ñ50), 
# and also in Yuan & Gomer (BJMSP 2020);
# ***--------------------------------------------------------------------***;
#====================== Table of Contents =======================#
# 1. Functions
# 2. Load data and run the functions
#================================================================#

######================== 1. Functions  =====================######                        

##==========================================##
# computing the sample skewness and 
# kurtosis with unidimensional data
##==========================================##                       
sk = function(x){
  n=length(x)
  x_bar=sum(x)/n
  x_c = scale(x, center = T, scale = F)
  x_c2=x_c^2
  m_x2 = sum(x_c2)/n # second-order moment
  m_x3 = sum(x_c2*x_c)/n # third-order moment
  m_x4 = sum(x_c2*x_c2)/n
  skew=m_x3/(m_x2*sqrt(m_x2))
  kurt=m_x4/(m_x2*m_x2)
  return(list(skew = skew, kurt = kurt))
} 

##==========================================##
# computing Mardia's multivariate skewness and kurtosis
##==========================================## 
Mardia = function(x){
  N = nrow(x) ; p = ncol(x)
  x_t = t(x)
  x_bar = colMeans(x)
  Scov = cov(x)*(N-1)/N
  scov_in = solve(Scov)
  mkurt = 0
    for(i in 1:N){
      x_i = x[i,]
      x_i0 = x_i - x_bar
      x_it0 = t(x_i0)
      d_i2 = x_it0%*%scov_in%*%x_i0
      mkurt = mkurt + d_i2*d_i2
    }
  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
   for(i in 1:N){
      x_i = x[i,]
      x_i0 = x_i - x_bar
      x_it0 = t(x_i0)
      d_i2 = x_it0%*%scov_in%*%x_i0
      mkurt = mkurt + d_i2*d_i2
      for(j in 1: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
      }
   }
  mskew = mskew/(N^2)
  
  mskew0 = N*mskew/6 #standardized Mardia's skewness
  df = p*(p+1)*(p+2)/6 # nominal degrees of freedom for standardized Mardia's skewness
  return(list(mskew0 = mskew0, mkurt0 = mkurt0, mkurt_p = mkurt_p, df = df))
  }

##==========================================##
# Robust transformation with Huber's type of weight
##==========================================##
robmusig = function(x, varphi=.05){
  ep = .00001 ; N = nrow(x) ; p = ncol(x);
  prob=1-varphi
  chip2 = qchisq(p = prob, df = p)
  kappa = sqrt(chip2)
  tau = (p*pchisq(chip2, df = p+2) + chip2*(1-prob))/p
  xr = matrix(0, nrow = N, ncol = p)
  
  # initial values
  sum_xt = rep(0, p) ; sum_xx = matrix(0, ncol = p, nrow = p)
  
  for(i in 1:N){
    x_it = t(x[i,])
    x_i = t(x_it)
    sum_xt = sum_xt + x_it
    sum_xx = sum_xx + x_i%*%as.matrix(x_it)
  }
  mu_t0 = as.matrix(sum_xt/N)
  sigma_0 = (sum_xx - N*t(mu_t0)%*%mu_t0)/(N - 1)
  
  delta = .01
  max.it = 500
  n_it = 1
  
  while (delta>ep && n_it <= max.it){
    sig_in = solve(sigma_0)
    sum_w1 = 0
    mu_t = rep(0, p)
    sigma = matrix(0, nrow = p, ncol = p)
    for(i in 1:N){
      x_it = t(x[i,])
      x_it0 = x_it - mu_t0
      x_i0 = t(x_it0)
      sigxi_0 = sig_in%*%x_i0
      d_i2 = x_it0%*%sigxi_0
      d_i = as.numeric(sqrt(d_i2))
      # Huber weight functions
      if(d_i<=kappa){
        w_i1 = 1
      }
      if(d_i>kappa){
        w_i1 = kappa/d_i
      }
      w_i2 = w_i1^2
      sum_w1 = sum_w1 + w_i1
      mu_t = mu_t + w_i1*x_it
      sigma = sigma + w_i2*x_i0%*%x_it0
    }
    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
    n_it = n_it + 1
  }
  if(n_it<max.it){
    err = 0
    hmu_t = mu_t1
    hsigma = sigma_1
    dw_vec = matrix(0, nrow = N, ncol = 4)
  for(i in 1:N){
    x_it = t(x[i,])
    x_it0 = x_it - mu_t0
    x_i0 = t(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){
      w_i1 = 1
    }
    if(d_i>kappa){
      w_i1 = kappa/d_i
    }
    w_i2 = w_i1%*%w_i1/tau
    xr[i,] = sqrt(w_i2)%*%x_it0
    dw_vec[i,] = c(i, d_i, w_i1, w_i2)
    }
  mean_xr = colMeans(xr)
  xr = xr + matrix(1, nrow = N, ncol = 1)%*%(hmu_t-mean_xr)
  dw_vec = dw_vec[order(dw_vec[,2], decreasing = T),] # order the case according to M-distance from large to small
 colnames(dw_vec) = c('i', 'd_i', 'w_i1', 'w_i2')
  id = t(1:N)
   }
  if (n_it>=max.it){
    err = 1
  print("The IRLS method of Huber-M does not converge within 500 iterations")
  }
  return(list(Robust_Mean = hmu_t, Robust_Covariance = hsigma, xr = xr, Error = err, dw_vec = dw_vec))
}

##==========================================##
# The following is the main program
##==========================================##
main = function(xmat, varphi = .05){
  # varphi: tuning parameter subject to change
  n = nrow(xmat) ; p = ncol(xmat)
  x = xmat
  # marginal skewness and kurtosis
  skew = rep(0, p)
  kurt = rep(0, p)
  skew_r = skew ; kurt_r = kurt
  for(j in 1:p){
    x_j = x[,j]
    temp = sk(x_j)
      skew[j] = temp$skew
      kurt[j] = temp$kurt
  }
  names(skew) = names(kurt) = 1:p
  #print("univariate skewness and kurtosis of raw data")
  #print(round(rbind(skew, kurt),3))
  
  Mardia.out = Mardia(as.matrix(x))
  #print("standardized Mardia's measures of skewness and kurtosis with raw data")
  #print(round(unlist(Mardia.out),3))
  
  out = list(univariate_skew.kurt_raw = round(rbind(skew, kurt),3), mardia_skew.kurt_raw = round(unlist(Mardia.out),3))
  temp = robmusig(varphi = varphi, x = as.matrix(x))
  #print(paste("Tuning parameter varphi with Huber-weights =", varphi))

  if(temp$Error==0){
    #print('Robust Means mu.hat')
    #print(round(temp$Robust_Mean,3))
    #print('Robust Covariance Matrix Sigma.hat')
    #print(round(temp$Robust_Covariance,3))

    for(j in 1:p){
      xr_j = temp$xr[,j]
      temp1 = sk(xr_j)
      skew_r[j] = temp1$skew # univariate skewness and kurtosis of robustified data
      kurt_r[j] = temp1$kurt
    }
    names(skew_r) = names(kurt_r) = 1:p
      #print("univariate skewness and kurtosis of robustified data")
      #print(round(rbind(skew_r, kurt_r),3))
  
    # standardized Mardia's measures of skewness and kurtosis with robustified data
    temp.Mardia = Mardia(x = temp$xr)
    Mardia.r.out = list(mskew_r = temp.Mardia$mskew0, mkurt_r = temp.Mardia$mkurt0, mkurt_rp = temp.Mardia$mkurt_p, df = temp.Mardia$df)
      #print("standardized Mardia's measures of skewness and kurtosis with robustified data")
      #print(round(unlist(Mardia.r.out),3))
  #print('Robustified Data')
  #print(temp$xr)
 out = list(univariate_skew.kurt_raw = round(rbind(skew, kurt),3), mardia_skew.kurt_raw = round(unlist(Mardia.out),3),
       univariate_skew.kurt_robust_data = round(rbind(skew_r, kurt_r),3), standardized_Mardia_skew.kurt_robust_data = round(unlist(Mardia.r.out),3), 
       Robust_Mean = round(temp$Robust_Mean,3), Robust_Covariance = round(temp$Robust_Covariance,3), Robustified_data = temp$xr)
  }
  return(out)
}

######================== 2. Load data and run the functions =====================######     
setwd() # set working directory
xmat = read.table("Mardia.dat") # read in data
main(xmat = xmat, varphi = .05) # run the main function