k <- 4
M <- 10
gamma <- matrix(rexp(k*M), nrow = M, ncol = k)
alpha <- rexp(k)
g_alpha <- c()
for (i in 1:k) {
g_alpha[i] <- M*(digamma(sum(alpha))-digamma(alpha[i]))
+ colSums(digamma(gamma)-digamma(rowSums(gamma)))[i]
}
diag(1,nrow = k)