r - 混合建模 - 来自 exp() 和 log() 的无限值的麻烦

Question

我正在为具有球面协方差结构的高斯混合模型编写一个函数——即$\Sigma_k = \sigma_k^2 I$. 此特定功能类似于mclust具有标识符 VII 的包。

http://en.wikipedia.org/wiki/Mixture_model

无论如何，我遇到的问题是权重矩阵的值无穷大。定义：令W是一个nxm矩阵，其中 n = 1, ..., n (obs 的数量) 和 m = 1, ..., m (mixtues 的数量)。W 的每个元素（即 w_ij）本质上可以定义为以下的特定形式：

w_im = \frac{a / b * exp(c)}{\sum_i=1^m [a_i / b_i * exp(c_i)]}

用数字计算这个给了我无限的价值。所以我正在尝试使用 log-identity log(x+y) = log(x) + log(1 + y/x)。但问题在于，它并不像log(x+y)而是log(\sum_i=1^m [a_i / b_i * exp(c_i)]).

这是一些代码定义：

1. n_im = a / b * exp(c);
2. d_.m = \sum_i=1^m [a_i / b_i * exp(c_i)]; 和

3. c_mat[i,j]作为第 [i,j] 项的指数值。

   n_mat[, i] <- log(a[i]) - log(b[i]) - c[,i] # numerator of w_im
   internal_vec1[i] <- (a[i] * b[1])/ (a[1] * b[i]) # an internal for the step below
   c_mat2 <- cbind(rep(1, n), c_mat[,1] - c_mat[,-1]) # since e^a / e^b = e^(a-b)
   for (i in 1:n) {
         d_vec[i] <- n_mat[i,1] + log(sum(internal_vec1 * exp(c_mat2[i,))) 
   } ## still getting infinite values
   

我试图尽可能简短地定义问题。整个函数显然比这大得多。但是，由于我遇到的问题是专门处理无限（和 1/infinity）值，我希望这个片段就足够了。有人在这里有编码技巧吗？

Answer 1

这是解决方案！（我在这上面花了太长时间了）

**第一个函数log_plus()解决了你想要 log(\sum_{i=1)^n x_i) 的简单问题

**第二个函数log_plus2()解决了上面描述的更复杂的问题，你想要 log(\sum_{i=1}^n [a_i / b_i * exp(c_i)])

log_plus <- function(xvec) {
  m <- length(xvec)
  x <- log(xvec[1])
  for (j in 2:m) {
    sum_j <- sum(xvec[1:j-1])
    x <- x + log(1 + xvec[j]/sum_j)
  }
  return(x)
}

log_plus2 <- function(a, b, c) {
  # assumes intended input of form sum(a/b * e^c)
  if ((length(a) != length(b)) || (length(a) != length(c))) {
    stop("Input equal length vectors")
  }
  if (!(all(c > 0) || all(c < 0))) {
    stop("All values of c must be either > 0 or  < 0.")
  }
  m <- length(a)
  # initilialize log sum
  x <- log(a[1]) - log(b[1]) + c[1]
  # aggregate / loop log sum
  for (j in 2:m) { 
    # build denominator
    b2 <- b[1:j-1]
    for (i in 1:j-1) {
      d1 <- 0
      c2 <- c[1:i]
      if (all(c2 > 0)) {
        c_min <- min(c2[1:j-1])
        c2 <- c2 - c_min 
      } else if (all(c2 < 0)) {
        c_min <- max(c2[1:j-1])
        c2 <- c2 - c_min
      }
      d1 <- d1 + a[i] * prod(b2[-i]) * exp(c2[i])
    }
    den <- b[j] * (d1)
    num <- a[j] * prod(b[1:j-1]) * exp(c[j] - c_min)
    x <- x + log(1 + num / den)
  }
  return(x)
}