2

我在使用 R 计算非平滑函数的梯度时遇到问题。

为了解释我的问题的核心，我创建了一个简化的问题。我有一个对数似然三变量正态密度函数，其中均值固定为零，方差的对角线项全为 1。

  library(mvtnorm)
  llk <- function(as,dat){

      sigma <-matrix(1,nrow=3,3)
      sigma[2,1] <- sigma[1,2] <-  as[1]
      sigma[3,1] <- sigma[1,3] <- as[2]
      sigma[3,2] <-  sigma[2,3] <- as[3]
      if (sum(eigen(sigma)$value < 0)>0) return(log(0))
      return(dmvnorm(dat,mean=c(0,0,0),sigma=sigma,log=TRUE))
      }

为了保持为 as 的所有值定义我的密度，如果方差不是正定矩阵，则 llk 会吐出 -Inf。现在，我想评估 llk 相对于 (0.5,0.2,0.3) 和 (0.9,-0.146,0.3) 处方差的三个非对角项的梯度。

使用 numDeriv 库中的 grad 函数，我可以将 (0.5,0.2,0.3) 处的梯度计算为：

  library(numDeriv)
  dat <- c(1,-1,2)
  as <- c(0.5,0.2,0.3)
  grad(llk,x=as,dat=dat)

 [1] -4.218858  4.533953 -6.128893

但是，此方法不适用于 (0.9,-0.146,0.3) 处的梯度。要看到这个：

  as <- c(0.9,-0.146,0.3)
  grad(llk,x=as,dat=dat)

 > Error in grad.default(llk, x = as, dat = dat) : 
   function returns NA at 9e-051.46e-053e-05 distance from x.

我知道这是因为第二个方差矩阵接近非正定。

有没有简单的方法来解决这个问题？谢谢你。

4

0 回答 0