R/One_way_median_polish.R

In hdftsa: High-Dimensional Functional Time Series Analysis

One_way_median_polish <-
function(Y, n_prefectures=51, year=1959:2020, age=0:100)
{
    #Number of years considered in the population
    n_year = length(year)
    #Number of ages considered in each year
    n_age = length(age)

    #######
    #Get the indexes for the row partitions by prefectures
    #######
    part_list = list()
    for(ik in 1:n_prefectures) {
      part_list[[ik]] = (n_year*ik-(n_year-1)):(n_year*ik)
    }
    row_partition_index=part_list


    if(length(class(Y)) == 2 && class(Y)[1] == "matrix" && class(Y)[2] == "array") 
    {
      sample_size <- n_year
      grid_size <- ncol(Y)
    } 
    else if(all(!inherits(Y, "matrix"), !inherits(Y, "array")))
    {
      stop("Y should be functional data")
    }

    Functional_grand_effect <- rep(0, length.out = grid_size)
    Functional_row_effect <- matrix(0, nrow = length(row_partition_index), ncol = grid_size)
    new_Functional_row_effect <- Functional_row_effect + 1
    Y_iter<-Y
    Iteration <- 1
    while (sum(new_Functional_row_effect) != 0) {
      message("Iteration=", Iteration,"\n")
      message("Sum =", sum(new_Functional_row_effect),"\n")

      partition_Y <- lapply(1:length(row_partition_index), function(k) {
        Y_iter[row_partition_index[[k]], ]})

      # Compute the median in each row with the curves in the row. Depth based selection

      Median_at_each_row<-lapply(partition_Y,get_median)

      # Get the functional median of the row functional medians
      Medians_medians_row<-get_median(list.rbind( Median_at_each_row))
      Functional_grand_effect=Functional_grand_effect+Medians_medians_row

      # Get the functional row effect
      new_Functional_row_effect<-t(sapply(Median_at_each_row, function(k) {k -Medians_medians_row}))
      # sum(new_Functional_row_effect)
      Functional_row_effect=Functional_row_effect+new_Functional_row_effect

      # Substract the row-median in each row from the rest of the curves in the row

      partition_Y<-lapply(partition_Y,Substraction)
      Y_iter<-list.rbind(partition_Y)
      Iteration <- Iteration + 1
    }
  return(list(grand_effect = Functional_grand_effect,
              row_effect = Functional_row_effect))
}