CODEWITHSUNDEEP
rmarkov-chainsmarkov
user2905427
user2905427

Reputation: 59

Use the markovchain package to compare two empirically estimated Markov chains

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)
}

Upvotes: 1

Views: 1393

Answers (1)

Giorgio Spedicato
Giorgio Spedicato

Reputation: 2503

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))

Upvotes: 3

Related Questions