coda gelman.diag():“ chol.default(W)中的错误:nn顺序的前导未成年人不是肯定的”

指法

当计算多元潜在比例缩减因子(MPSRF)时gelman.diag()coda包中函数将引发错误。

> load("short_mcmc_list.rda")
> niter(short.mcmc.list)
[1] 100
> nvar(short.mcmc.list)
[1] 200
> nchain(short.mcmc.list)
[1] 2
> 
> coda::gelman.diag(
    short.mcmc.list, 
    autoburnin = FALSE,
    multivariate = TRUE
)

chol.default(W)中的错误:199阶的前导未成年人不是肯定的

这个错误是什么意思?

此问题先前发布在R coda上,“第3阶的未成年人不是肯定的”主要结论是:“结论:获得Gelman-Rubin诊断的多元估计值似乎有问题。设置多元= FALSE可以解决问题,并为每个变量输出单个变量估计。” 它已经6岁了,所以答案可能已过时。

指法

在中gelman.diag(),MPSRF的计算公式为:

> coda::gelman.diag <-
function (x, confidence = 0.95, transform = FALSE, autoburnin = FALSE,
        multivariate = TRUE)
{
#A lot of code ...

Niter <- niter(x)
Nchain <- nchain(x)
Nvar <- nvar(x)
#W is the variance matrix of the chain
S2 <- array(sapply(x, var, simplify = TRUE), dim = c(Nvar,
                                                   Nvar, Nchain))
W <- apply(S2, c(1, 2), mean)
#A lot of code ...

if (Nvar > 1 && multivariate) {
CW <- chol(W)
#This is max eigenvalue from W^-1*B.
emax <- eigen(
  backsolve(CW, t(backsolve(CW, B, transpose = TRUE)), transpose = TRUE),
  symmetric = TRUE, only.values = TRUE)$values[1]
}

gelman.diag()删除了产生错误的Cholesky分解,并将W和B添加到要返回的列表中,从而对进行了编辑。这使我明白了为什么W不能进行Cholesky分解。

my.gelman.diag <- function(x,
                           confidence = 0.95,
                           transform = FALSE, 
                           autoburnin = FALSE,
                           multivariate = TRUE
){
  x <- as.mcmc.list(x)
  if (nchain(x) < 2)
    stop("You need at least two chains")
  if (autoburnin && start(x) < end(x)/2)
    x <- window(x, start = end(x)/2 + 1)
  Niter <- niter(x)
  Nchain <- nchain(x)
  Nvar <- nvar(x)
  xnames <- varnames(x)
  if (transform)
    x <- gelman.transform(x)
  x <- lapply(x, as.matrix)
  S2 <- array(sapply(x, var, simplify = TRUE),
              dim = c(Nvar, Nvar, Nchain)
  )
  W <- apply(S2, c(1, 2), mean)
  xbar <- matrix(sapply(x, apply, 2, mean, simplify = TRUE),
                 nrow = Nvar, ncol = Nchain)
  B <- Niter * var(t(xbar))
  if (Nvar > 1 && multivariate) {  #ph-edits 
    # CW <- chol(W)
    #    #This is W^-1*B.
    # emax <- eigen(
    #  backsolve(CW, t(backsolve(CW, B, transpose = TRUE)), transpose = TRUE),
    # symmetric = TRUE, only.values = TRUE)$values[1]
    emax <- 1
    mpsrf <- sqrt((1 - 1/Niter) + (1 + 1/Nvar) * emax/Niter)
  }  else {
    mpsrf <- NULL
  }

  w <- diag(W)
  b <- diag(B)
  s2 <- matrix(apply(S2, 3, diag), nrow = Nvar, ncol = Nchain)
  muhat <- apply(xbar, 1, mean)
  var.w <- apply(s2, 1, var)/Nchain
  var.b <- (2 * b^2)/(Nchain - 1)
  cov.wb <- (Niter/Nchain) * diag(var(t(s2), t(xbar^2)) - 2 *
                                    muhat * var(t(s2), t(xbar)))
  V <- (Niter - 1) * w/Niter + (1 + 1/Nchain) * b/Niter
  var.V <- ((Niter - 1)^2 * var.w + (1 + 1/Nchain)^2 * var.b +
              2 * (Niter - 1) * (1 + 1/Nchain) * cov.wb)/Niter^2
  df.V <- (2 * V^2)/var.V
  df.adj <- (df.V + 3)/(df.V + 1)
  B.df <- Nchain - 1
  W.df <- (2 * w^2)/var.w
  R2.fixed <- (Niter - 1)/Niter
  R2.random <- (1 + 1/Nchain) * (1/Niter) * (b/w)
  R2.estimate <- R2.fixed + R2.random
  R2.upper <- R2.fixed + qf((1 + confidence)/2, B.df, W.df) *
    R2.random
  psrf <- cbind(sqrt(df.adj * R2.estimate), sqrt(df.adj * R2.upper))
  dimnames(psrf) <- list(xnames, c("Point est.", "Upper C.I."))
  out <- list(psrf = psrf, mpsrf = mpsrf, B = B, W = W) #added ph
  class(out) <- "gelman.diag"
  return( out )
}

适用my.gelman.diag()short.mcmc.list

> l <- my.gelman.diag(short.mcmc.list, autoburnin = FALSE, multivariate = TRUE)
> W <- l$W  #Within-sequence variance
> B <- l$B  #Between-sequence variance

对W的研究表明W确实是正定的,但其特征值接近0,因此失败。

> evals.W <- eigen(W, only.values = TRUE)$values
> min(evals.W)
[1] 1.980596e-16

确实,增加容差表明W确实是正定的。

> matrixNormal::is.positive.definite(W, tol = 1e-18)
[1] TRUE

所以实际上,W接近线性相关性...

> eval <- eigen(solve(W)%*%B, only.values = TRUE)$values[1]

resolve.default(W)中的错误:系统在计算上是单数的:倒数条件数= 7.10718e-21

因此,实际上,删除W的最后两列使其成为线性独立/正定的。这表明链中存在相关的参数,并且可以减少参数的数量。

> evals.W[198]
[1] 9.579362e-05
> matrixNormal::is.positive.definite(W[1:198, 1:198])
[1] TRUE
> chol.W <- chol(W)

W是马尔可夫链的内方差,衡量状态中每个值与均值的差异。如果W接近奇异,则方差很小,因此链条变化不大。这是一条缓慢移动的链条。用户应使用轨迹图进行调查,并可能减少参数数量。链条也可能太短;如果链较长,则链中的值可能会足够不同,以使W线性独立。

为了避免函数崩溃,我建议使用purrr::possibly()替换缺少的值,而不是引发过时的错误。

>  pgd <- purrr::possibly(coda::gelman.diag, list(mpsrf = NA), quiet = FALSE)
> pgd(short.mcmc.list, autoburnin = FALSE, multivariate = TRUE)
Error: the leading minor of order 199 is not positive definite
[1] NA

有关更多信息,请参见Brooks和Gelman,1992年。

参考:Gelman,A和Rubin,DB(1992)使用多个序列进行迭代模拟的推论,《统计科学》,第7卷,第457-511页。

本文收集自互联网,转载请注明来源。

如有侵权,请联系 [email protected] 删除。

编辑于
0

我来说两句

0 条评论
登录 后参与评论

相关文章

chol.default(Cxx)中的错误:顺序的前导未成年人不是肯定的

在chol.default(K)中出现错误:5级前导未成年人不是使用betareg的正定

R 中 MCMCglmm 模型的 Gelman-Rubin 统计量

如何更改未成年人的数量?

R 中 coda 中 mcmc 函数的细化

.hbs语法在Coda(OSX)中的突出显示

在MATLAB中替代diag(X'* C * X)

iBeacon约束可以接受多个专业/未成年人

我怎么知道Linux模块初始化的未成年人

SQL-未成年人之前存在主要帐户

如何更改Coda 2中的全局字体大小?

我正在尝试在python中创建一个程序,以接受用户提供的输入并告诉他们是否未成年人

R中的chol()函数不断返回上三角(我需要下三角)

np.diag 不包括数组中的 0

尝试在春季启动时集成Coda Hale指标,而不是/ metrics中的指标

如何在Spring Boot应用程序中配置HikariCP和Dropwizard / Coda-Hale指标

ArrayFire.jl 和 chol

plot_model(type =“ diag”)中的错误:参数隐含不同的行数错误

厨师距离的dotplot_diag(lmer模型)-主题值而不是索引

为什么diag功能这么慢?[在R 3.2.0或更早版本中]

在OpenMDAO的ExecComp中,shape_by_conn与has_diag_partials兼容吗?

如何在tensorflow1.x中删除二维张量的下三角(包括diag)?

numpy diag不返回矩阵

R chol2inv() 方法给了我奇怪的结果

当在hmmlearn包中使用model.score时,为什么会出现“'diag'混合Covar必须为非负数”错误?

DIAG [S1000] [SAP AG][LIBODBCHDB SO][HDBODBC] 一般错误;-10427 参数/列 (8) 从数据类型 NVARCHAR 到 ASCII 的转换失败

为什么Numpy Diag函数的行为很奇怪?

1./diag(A) 是什么意思?

紧凑/高效地替代diag(XVX ^ T)?