使用 markovchain 包比较两个经验估计的马尔可夫链

Use the markovchain package to compare two empirically estimated Markov chains

我需要比较两个概率矩阵以了解链的接近程度,因此我将使用测试的结果 P 值。

我尝试使用 markovchain r 包,更具体地说是 divergenceTest 函数。但是,问题是该功能没有正确实现。它是基于第139页“Statistical Inference Based on Divergence Measures”一书的测试,我联系了包开发人员,但他们仍然没有纠正,所以我试图实现,但我遇到了麻烦,谁能帮助我找到错误?

参数:freq_matrix:是一个频率矩阵,用来估计概率矩阵。 hypothetic: 是用来和估计矩阵比较的矩阵。

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

发散测试已被 verifyHomogeneity 功能取代。它需要并输入可以强制转换为原始转换矩阵的元素列表(从 createSequenceMatrix 开始)。然后测试它们是否属于同一个未知的 DTMC。

参见下面的示例:

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