我不知道如何使用markovchain包来实现这一点，但您可以编写一个函数来产生类似的结果，下面的函数通过MLE生成转移概率的估计值，但标准误差置信区间的上限和下限。你可以使用标准误差和估计值来建立正态理论置信区间（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)

1.创建函数来估计模型并打印结果。

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

1.使用我生成的函数估计模型：

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++编写的，我还没有研究源代码。如果这个差异是不可接受的，也许其他人能更接近。