在R中，如何手动编写正态分布核的似然性？

frebpwbc 于 2023-02-20 发布在其他

关注(0)|答案(1)|浏览(131)

具体来说，x和mu之差的乘积，精度矩阵，x和mu之差的转置，你是怎么编码的，我下面的代码正确吗？提前感谢。

(colSums(dat-mu_mat)%*%solve(sigma)%*%colSums(dat-mu_mat))

其中mu_mat是向下重复n次的均值的行向量。
完整代码如下：

dat<-rmvnorm(100,mean=c(200,0.1),sigma=matrix(c(5,0,0,0.02),nrow=2))
n<-nrow(dat)
mu<-matrix(c(200,0.1),nrow=1)
mu_mat<-matrix(rep(c(200,0.1),100),nrow=100,ncol=2,byrow=TRUE)

loglik_mvn<-function(n,d,x,mu_mat,sigma){
  (-n*d/2)*log(2*pi)-(n/2)*det(sigma,log=TRUE)-0.5*(colSums(x-mu_mat)%*%solve(sigma)%*%colSums(x-mu_mat))

loglik_mvn(nrow(dat),ncol(dat),dat,mu_mat,sigma)
dmvnorm(dat,mean=mu,sigma=sigma,log=TRUE)
}

r

来源：https://stackoverflow.com/questions/75479505/how-do-you-code-the-likelihood-of-the-normal-distribution-kernel-in-r-manually

1条答案

按热度按时间

qoefvg9y1#

如果您有数据集、均值向量和sigma，则可以计算指数部分：

sum(diag(solve(sig, tcrossprod(t(x) - mu))))

因此，您可能会看起来像：

log_like <- function(x, mu, sig){
  -sum(diag(solve(sig, tcrossprod(t(x) - mu))))/2 - #The part you are interested in
    nrow(x)*(ncol(x)*log(2*pi) + log(det(sig)))/2
}

x <- iris[-5]
mu <- colMeans(x)
sig <- var(x)

log_like(x, mu, sig)
#> [1] -379.9213

sum(mvtnorm::dmvnorm(x, mu, sig, TRUE)) # Produces same results as above
#> [1] -379.9213

创建于2023年2月16日，使用reprex v2.0.2
您还可以使用以下命令：

log_like2 <- function(x, mu, sig){
  u <- t(x) - mu #Notice mu is a vector and not a matrix
  -sum(u*solve(sig, u))/2 -  # the part you are interested in
   nrow(x)*(ncol(x)*log(2*pi) + log(det(sig)))/2
}

log_like2(x, mu, sig)
#> [1] -379.9213

把所有的东西放在一起，利用R的循环：

log_like3 <- function(x, mu, sig){
  u <- t(x) - mu
  -sum(u * solve(sig, u) + log(2*pi) + log(det(sig))/ncol(x))/2
}
log_like3(x, mu, sig)
[1] -379.9213

赞(0）回复(0）举报 2023-02-20

我来回答

在R中，如何手动编写正态分布核的似然性？

1条答案

相关问题

热门标签

最新问答