我在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) )
。
1条答案
按热度按时间cnjp1d6j1#
在您代码中
与更紧凑的方法相比:
虽然这两个表达式给出了相似的结果,但第一个表达式似乎比第二个表达式更快。我不知道为什么,因为第一个表达式需要计算n * n矩阵。for循环需要for-ever来计算。
一个二个一个一个
NB:我花了一段时间才得到您使用的种子。下次将其包含在您的数据生成中