I am using HSMM(Hidden Semi Markov Model) R code with my fMRI data files. However, if I reverse the order of the files and run the code, the results(state correlation matrix, state transition probability etc.) keep changing because of the initialization with E-M algorithm issue. I want to get consistent results even if the order of the files is reversed, is there a solution?

library(depmixS4)
library(lattice)
library(MASS)
library(remotes)
library(distr)
library(mvnmle)
library(mvtnorm)
library(mhsmm)
library(gplots)
library(rasterVis)
library(lattice) 
library(shapes)

setwd("/Users/jarvis/Desktop/MyDataSet/")

filenames <- Sys.glob(file.path(getwd(),'*','sub-*','sub-*Average_ROI_n50.csv'))

ilenames

DataReadIn <- lapply(filenames,function(i){
  table<-read.csv(i, header = TRUE, row.names=NULL)
  table<-scale(table)
  return(table)
})

DataFinal<-do.call(rbind, DataReadIn)
DataFinal<-as.data.frame(DataFinal)

TimeSerLen <- lapply(filenames,function(i){
  table<-read.csv(i, header = TRUE, row.names=NULL)
  Nrows<-nrow(table)
  return(Nrows)
})

Lengths<-do.call(rbind, TimeSerLen)

number_states = 7

number_runs <- nrow(Lengths)
number_regions <- ncol(DataFinal)

output_array = list()
seglength = ceiling(min(Lengths[, 1]) / number_states)

for (i in 1:number_states) 
{
  output_array[[i]] <- matrix(NA, nrow = 1, ncol = number_regions)
  colnames(output_array[[i]]) <- colnames(DataFinal)
  for(j in 1:number_runs){
    output_array[[i]] <- rbind(output_array[[i]], DataFinal[(1+sum(Lengths[1:(j - 1), 1]) + seglength * (i - 1)):(sum(Lengths[1:(j - 1), 1]) + seglength * i), ])
  }
  output_array[[i]] <- output_array[[i]][-1, ]
}

MeanEst <- lapply(output_array,function(i){
  means <- mlest(i)$muhat
  return(means)
})

CovEst <- lapply(output_array,function(i){
  covs <- mlest(i)$sigmahat
  return(covs)
})

initial <- rep(1/number_states, number_states)

TransMatrix <- matrix(NA, nrow=number_states, ncol=number_states)
diag(TransMatrix) <- 0
TransMatrix[is.na(TransMatrix)] <- (1/(number_states-1))

b <- list(mu = MeanEst, sigma = CovEst)

Datanew <- list(x = DataFinal, N = Lengths[,1])
class(Datanew) <- "hmm.data"

d <- list(shape = sample(1:7, size=number_states), scale = sample(10:50, size=number_states), type = "gamma")

model1 <- hsmmspec(init = initial, trans = TransMatrix, parms.emis = b, dens.emis = dmvnorm.hsmm, sojourn = d)

testrun <- hsmmfit(Datanew, model1,mstep=mstep.mvnorm, lock.transition=FALSE, maxit=1000, lock.d=FALSE)