Use the markovchain package to compare two empirically estimated Markov chains

Question

I need to compare two probability matrices to know the degree of proximity of the chains, so I would use the resulting P-Value of the test.

I tried to use the markovchain r package, more specifically the divergenceTest function. But, the problem is that the function is not properly implemented. It is based on the test of the book "Statistical Inference Based on Divergence Measures" on page 139, I contacted the package developers, but they still have not corrected, so I tried to implement, but I'm having trouble, could anyone help me to find the error?

Parameters: freq_matrix: Is a frequency matrix used to estimate the probability matrix. hypothetic: Is the matrix used to compare with the estimated matrix.

divergenceTest3 <- function(freq_matrix, hypothetic){  
  n <- sum(freq_matrix)
  empirical = freq_matrix
  for (i in 1:length(hypothetic)){
    empirical[i,] <- freq_matrix[i,]/rowSums(freq_matrix)[i]
  }
  M <- nrow(empirical)
  v <- numeric()
  out <- 2 * n / .phi2(1)
  sum <- 0
  c <- 0  
  for(i in 1:M){    
    sum2 <- 0
    sum3 <- 0    
    for(j in 1:M){
      if(hypothetic[i, j] > 0){
        c <- c + 1
      }      
      sum2 <- sum2 + hypothetic[i, j] * .phi(empirical[i, j] / hypothetic[i, j])
    }    
    v[i] <- rowSums(freq_matrix)[i]
    sum <- sum + ((v[i] / n) * sum2)
  }
  TStat <- out * sum
  pvalue <- 1 - pchisq(TStat, c-M)  
  cat("The Divergence test statistic is: ", TStat, " the Chi-Square d.f. are: ", c-M," the p-value is: ", pvalue,"\n")
  out <- list(statistic = TStat, p.value = pvalue)  
  return(out)
}
# phi function for divergence test
.phi <- function(x) {
  out <- x*log(x) - x + 1
  return(out)
}
# another phi function for divergence test
.phi2 <- function(x) {
  out <- 1/x
  return(out)
}

score 3 · Accepted Answer · answered Feb 20 '17 at 00:20

The divergence test has been replaced by the verifyHomogeneityfunction. It requires and input list of elements that can be coerced to a raw transition matrix (as of createSequenceMatrix). Then it tests whether they belong to the same unknown DTMC.

See the example below:

myMatr1<-matrix(c(0.2,.8,.5,.5),byrow=TRUE, nrow=2)
myMatr2<-matrix(c(0.5,.5,.4,.6),byrow=TRUE, nrow=2)
mc1<-as(myMatr1,"markovchain")
mc2<-as(myMatr2,"markovchain")
mc
mc2
sample1<-rmarkovchain(n=100, object=mc1)
sample2<-rmarkovchain(n=200, object=mc2)
# should reject
verifyHomogeneity(inputList = list(sample1,sample2))
#should accept
sample2<-rmarkovchain(n=200, object=mc1)
verifyHomogeneity(inputList = list(sample1,sample2))

Hello professor, this function seems to work fine now! But, in the documentation only says that the test of article [S. Kullback at all, Tests for Contingency Tables and Markov Chains], correct? I would like to know specifically which equation you have implemented, since I need to understand the test and I am having difficulties. Thanks in advance! — user2905427, May 03 '17 at 07:52

Use the markovchain package to compare two empirically estimated Markov chains

1 Answers1