R语言 实现使用矩阵运算计算对数似然的特定方法

y3bcpkx1  于 2023-02-20  发布在  其他
关注(0)|答案(1)|浏览(166)

我在CrossValidated.com answer中遇到了对数似然的数学表达式,不清楚应该如何在R中实现。我不确定SO是否可以像CV一样表示MathML,但这是第二个(未接受)答案中的第一个等式:

$$
\begin{eqnarray}
    \ell(\mu, \Sigma) &=& C - \frac{m}{2}\log|\Sigma|-\frac{1}{2} \sum_{i=1}^m \text{tr}\left[(\mathbf{x}^{(i)}-\mu)^T \Sigma^{-1} (\mathbf{x}^{(i)}-\mu)\right]\\
$$

我关注的是该等式中的第三项,根据该页上的另一个答案,我认为跟踪操作是不必要的。我想我可以查看现有各种包中的几个实现之一,但我认为它们使用了更经济的方法,这些方法没有明确遵循该等式的过程,如answer here:中的@onyambu
我从前面的一个SO示例中删除代码:

library(MASS)
# Make covariance matrix. See note above re the implications of using a correlation matrix.
S = matrix(c(1.0, 0.2, 0.1, 0.35, 0.0,
0.2, 1.0, 0.0, 0.4, 0.0,
0.1, 0.0, 1.0, 0.0, 0.4,
0.35, 0.4, 0.0, 1.0, 0.6,
0.0, 0.0, 0.4, 0.6, 1.0), ncol = 5)
colnames(S) = c("Y1", "X1", "X2", "Z1" ,"Z2")
rownames(S) = colnames(S)
# Make mean vector
mus = c(1, 2, 3, 4, 5); names(mus) = colnames(S)
# Generate 5347 observations
obs = mvrnorm(n = 200, mu = mus, Sigma = S)

这是对一个问题的回答,但没有使用矩阵表达式的求和,我想我可以用一个for循环为每个数据点创建单独的贡献:

llmat.term3 <- matrix(NA, 200,1)
for(n in 1:200) {
     llmat.term3[n] <- t(obs[n,]-mus) %*% solve(S) %*% (obs[n,]-mus) }
sum(llmat.term3)
#[1] 982.7356

...但我想知道是否有一个更紧凑的矩阵方法?或者我想,填补了空白,我的线性代数知识,解释为什么sum(u * solve(sig, u)是相同的sum{i=1,N} ( t(obs[n,]-mu) %*% S^-1 %*% (obs[n,]-mu) )

cnjp1d6j

cnjp1d6j1#

在您代码中

S = matrix(c(1.0, 0.2, 0.1, 0.35, 0.0,
             0.2, 1.0, 0.0, 0.4, 0.0,
             0.1, 0.0, 1.0, 0.0, 0.4,
             0.35, 0.4, 0.0, 1.0, 0.6,
             0.0, 0.0, 0.4, 0.6, 1.0), ncol = 5)
colnames(S) = c("Y1", "X1", "X2", "Z1" ,"Z2")
rownames(S) = colnames(S)
# Make mean vector
mus = c(1, 2, 3, 4, 5); names(mus) = colnames(S)
# Generate 5347 observations
set.seed(123)
obs = MASS::mvrnorm(n = 200, mu = mus, Sigma = S)
llmat.term3 <- matrix(NA, 200,1)
for(n in 1:200) {
  llmat.term3[n] <- t(obs[n,]-mus) %*% solve(S) %*% (obs[n,]-mus) }
sum(llmat.term3)
#[1] 982.7356

与更紧凑的方法相比:

u <- t(obs) - mus

sum(diag(solve(S, tcrossprod(u))))
#> [1] 982.7356    
sum(u * solve(S, u))
#> [1] 982.7356

虽然这两个表达式给出了相似的结果,但第一个表达式似乎比第二个表达式更快。我不知道为什么,因为第一个表达式需要计算n * n矩阵。for循环需要for-ever来计算。
一个二个一个一个
NB:我花了一段时间才得到您使用的种子。下次将其包含在您的数据生成中

相关问题