使用 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))
我需要比较两个概率矩阵以了解链的接近程度,因此我将使用测试的结果 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))