logo标签列表

离散时间马尔可夫链的标准误差和标准差计算

rmatrixstandard-error
3
我有一个从一种状态转移到另一种状态的计数矩阵，我想要计算最大似然估计，标准误差和标准偏差。 "markovchain"软件包有一个示例，但数据是一个序列。我的数据来自于155家公司的平衡面板数据集，所以他们提供的示例代码对我无效。
这是我遵循的示例：
 data(rain)
 rain$rain[1:10]

 [1] "6+"  "1-5" "1-5" "1-5" "1-5" "1-5" "1-5" "6+"  "6+"  "6+"

 #obtaining the empirical transition matrix
 createSequenceMatrix(stringchar = rain$rain)
      0 1-5  6+ 
0   362 126  60 
1-5 136  90  68 
6+   50  79 124

#fitting the DTMC by MLE
alofiMcFitMle <- markovchainFit(data = rain$rain, method = "mle", name = "Alofi")
alofiMcFitMle
$estimate
Alofi
A  3 - dimensional discrete Markov Chain defined by the following states:
 0, 1-5, 6+
 The transition matrix  (by rows)  is defined as follows:
            0       1-5        6+
0   0.6605839 0.2299270 0.1094891
1-5 0.4625850 0.3061224 0.2312925
6+  0.1976285 0.3122530 0.4901186
$standardError
             0        1-5         6+
0   0.03471952 0.02048353 0.01413498
1-5 0.03966634 0.03226814 0.02804834
6+  0.02794888 0.03513120 0.04401395
$confidenceInterval
$confidenceInterval$confidenceLevel
[1] 0.95
$confidenceInterval$lowerEndpointMatrix
            0       1-5         6+
0   0.6034754 0.1962346 0.08623909
1-5 0.3973397 0.2530461 0.18515711
6+  0.1516566 0.2544673 0.41772208
$confidenceInterval$upperEndpointMatrix
            0       1-5        6+
0   0.7176925 0.2636194 0.1327390
1-5 0.5278304 0.3591988 0.2774279
6+  0.2436003 0.3700387 0.5625151
$logLikelihood
[1] -1040.419

因为我已经有了一个计数数据矩阵，所以无法使用上述代码。我只想采用我的6x6转换计数矩阵来确定最大似然估计值、标准误差（置信区间）和标准偏差。有没有人有一个我可以跟随的例子？

我没有看到一个合适的统计请求。你提供了一组单一的观察结果，显然是某个过程的结束状态，并且正在要求使用机器学习方法完成 <something>，这与未观察到的转换过程有关。 - IRTFM
是的。我想计算一个6x6矩阵中每个转移的MLE，即矩阵中每个元素的计数除以行总和。一旦我有了MLE，我想获得MLE的标准误差和标准偏差。 - gm007
仍不清楚MLE将是什么估计。MLE通常指统计模型的参数。如果您尚未观察到完整的过程，则无法计算出一个参数。我能想象从单个6x6矩阵中产生的最佳结果是比例的估计值，se可能来自多项式分布。然而，并不能保证这些估计值是转换矩阵的估计值，因为我们不知道起始点以及该过程是否已达到任何稳定状态。 - IRTFM
我应该提到我所遵循的过程。抱歉。我假设一个一阶马尔可夫链。我有6个状态，代表了样本期间155家公司的信用评级。我计算从状态i开始的公司数量，并计算转移到状态j的公司数量以制作矩阵。估计值是从状态i转移到状态j的公司比例，矩阵是系统的概率分布。 - gm007
我不明白为什么一个具有未确定起始状态的马尔可夫过程的观察计数必须提供足够的信息来估计过程转移矩阵。在我看来，您没有提供足够的信息。 - IRTFM
1个回答
1
我不确定如何使用markovchain包来完成此操作，但您可以编写一个函数来产生类似的结果。下面的函数通过MLE生成转换概率的估计值，但通过参数化bootstrap生成标准误差、下限和上限置信区间。您可以使用标准误差和估计值制作正常理论置信区间(ci = "normal")，也可以使用参数化bootstrap分布的相关分位数(ci = "quantile")。 parallel标志标识矩阵的级别是否应被视为有序。如果有序，则模型使用有序登录来计算概率。如果无序(即parallel = FALSE)，则使用多项式登录来计算概率。
如果我理解正确，您正在从createSequenceMatrix(stringchar = rain$rain)输出开始，而没有生成它的基础值序列。以下是如何操作的步骤:
  1. 生成数据。
library(markovchain)
#> Package:  markovchain
#> Version:  0.9.0
#> Date:     2022-07-01
#> BugReport: https://github.com/spedygiorgio/markovchain/issues
data(rain)
mat <-  createSequenceMatrix(stringchar = rain$rain)

制作估计模型并打印结果的功能。
mcmat_est <- function(mat, parallel = FALSE, ci = c("normal", "quantile"), conf_level=.95, R=2500){
  ci <- match.arg(ci)
  if(!is.null(rownames(mat))){
    vals <- rownames(mat)
  }else{
    vals <- 1:nrow(mat)
  }
  state_vals <- lapply(vals, \(i)data.frame(vals = rep(vals, mat[i, ])))
  if(parallel){
    state_vals <- lapply(state_vals, \(x){x$vals <- factor(x$vals, levels=vals, ordered=TRUE); x})
    mods <- lapply(state_vals, \(x)MASS::polr(vals ~ 1, data=x, Hess = TRUE))
    draws <- lapply(mods, \(m)cbind(-Inf, MASS::mvrnorm(R, m$zeta, vcov(m)), Inf))
    qs <- lapply(draws, \(d)plogis(d))
    probs<- lapply(qs, \(x)x[,-1] - x[,-ncol(x)])
    ests <- lapply(mods, \(m)plogis(c(-Inf, m$zeta, Inf)))
    ests <- lapply(ests, \(e)e[-1]-e[-length(e)])
  }else{
    state_vals <- lapply(state_vals, \(x){x$vals <- factor(x$vals, levels=vals); x})
    mods <- lapply(state_vals, \(x)nnet::multinom(vals ~ 1, data=x, Hess = TRUE, trace=FALSE))
    draws <- lapply(mods, \(m)MASS::mvrnorm(R, c(0, coef(m)), 
                                            cbind(0, rbind(0, vcov(m)))))
    probs <- lapply(draws, \(d)prop.table(exp(d), 1))
    ests <- lapply(mods, \(m)exp(c(0, coef(m))))
    ests <- lapply(ests, \(e)e/sum(e))
  }
  logLik <- sum(sapply(mods, logLik))
  nobs <- sum(sapply(mods, \(x)attr(logLik(x), "nobs")))
  df <- sum(sapply(mods, \(x)attr(logLik(x), "df")))
  attr(logLik, "nobs") <- nobs
  attr(logLik, "df") <- df
  class(logLik) <- "logLik"
  est <- do.call(rbind, ests)
  se <- do.call(rbind, lapply(probs, \(p)apply(p, 2, sd)))
  crit_z <- 1-(1-conf_level)/2
  if(ci == "normal"){
    lwr <- est - qnorm(crit_z)*se
    upr <- est + qnorm(crit_z)*se
  }else{
    lwr <- do.call(rbind, lapply(probs, \(p)apply(p, 2, quantile, 1-crit_z)))
    upr <- do.call(rbind, lapply(probs, \(p)apply(p, 2, quantile, crit_z)))
  }
  rownames(est) <- rownames(se) <- rownames(lwr) <- rownames(upr) <- rownames(mat)
  colnames(est) <- colnames(se) <- colnames(lwr) <- colnames(upr) <- colnames(mat)
  res <- list(estimate = est, se = se, lower=lwr, upper=upr, logLik=logLik)
  class(res) <- "mcest"
  res
}
print.mcest <- function(x, ..., digits=3){
  cat("\nEstimated Transition Probabilities\n")
  print(x$estimate, digits=digits)
  cat("\nStandard Errors of Transition Probabilities\n")
  print(x$se, digits=digits)
  cat("\nLower Confidence Bounds\n")
  print(x$lower, digits=digits)
  cat("\nUpper Confidence Bounds\n")
  print(x$upper, digits=digits)
  cat('\n')
  print(x$logLik)
}

使用我生成的函数估计模型：
out <- mcmat_est(mat)  
out
#> 
#> Estimated Transition Probabilities
#>         0   1-5    6+
#> 0   0.661 0.230 0.109
#> 1-5 0.463 0.306 0.231
#> 6+  0.198 0.312 0.490
#> 
#> Standard Errors of Transition Probabilities
#>          0    1-5     6+
#> 0   0.0203 0.0181 0.0135
#> 1-5 0.0281 0.0269 0.0248
#> 6+  0.0245 0.0289 0.0309
#> 
#> Lower Confidence Bounds
#>         0   1-5    6+
#> 0   0.621 0.194 0.083
#> 1-5 0.407 0.253 0.183
#> 6+  0.150 0.256 0.430
#> 
#> Upper Confidence Bounds
#>         0   1-5    6+
#> 0   0.700 0.265 0.136
#> 1-5 0.518 0.359 0.280
#> 6+  0.246 0.369 0.551
#> 
#> 'log Lik.' -1040.419 (df=6)
  1. 与 markovchainFit() 的输出进行比较。
fit <- markovchainFit(data = rain$rain, method = "mle", name = "Alofi")
fit
#> $estimate
#> Alofi 
#>  A  3 - dimensional discrete Markov Chain defined by the following states: 
#>  0, 1-5, 6+ 
#>  The transition matrix  (by rows)  is defined as follows: 
#>             0       1-5        6+
#> 0   0.6605839 0.2299270 0.1094891
#> 1-5 0.4625850 0.3061224 0.2312925
#> 6+  0.1976285 0.3122530 0.4901186
#> 
#> 
#> $standardError
#>              0        1-5         6+
#> 0   0.03471952 0.02048353 0.01413498
#> 1-5 0.03966634 0.03226814 0.02804834
#> 6+  0.02794888 0.03513120 0.04401395
#> 
#> $confidenceLevel
#> [1] 0.95
#> 
#> $lowerEndpointMatrix
#>             0       1-5        6+
#> 0   0.5925349 0.1897800 0.0817850
#> 1-5 0.3848404 0.2428780 0.1763188
#> 6+  0.1428496 0.2433971 0.4038528
#> 
#> $upperEndpointMatrix
#>             0       1-5        6+
#> 0   0.7286330 0.2700740 0.1371931
#> 1-5 0.5403296 0.3693669 0.2862663
#> 6+  0.2524073 0.3811089 0.5763843
#> 
#> $logLikelihood
#> [1] -1040.419

本文由reprex package (v2.0.1)创建于2022年12月15日

这两个模型的对数似然相同，这表明其基础估计过程相似。但标准误差不同 - 原始函数没有使用参数自助法。我不确定它们做了什么，因为用于估计模型的函数是用C++编写的，我还没有研究源代码。如果这种差异是不可接受的，也许其他人能够更接近。

@ DaveArmstrong：非常感谢您的回答！我需要仔细阅读并在有问题时与您联系！ - stats_noob
另外，我一直在研究这个问题 - 你能否请看一下并告诉我是否我已经正确地完成了这个问题？谢谢！https://stackoverflow.com/questions/74806562/correctly-using-clustereval-in-r - stats_noob
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接