我正在尝试在R中实现多元正态分布的一阶导数,基于这里和这里发布的多元正态分布的rcpp
实现。
下面是一个快速的R实现:
mvnormDeriv = function(..., mu=rep(0,length(list(...))), sigma=diag(length(list(...)))) {
if(sd(laply(list(...),length))!=0)
stop("The vectors not same length.")
fn = function(x) -1 * c((1/sqrt(det(2*pi*sigma))) * exp(-0.5*t(x-mu)%*%solve(sigma)%*%(x-mu))) * solve(sigma,(x-mu))
out = t(apply(cbind(...),1,fn))
colnames(out) = c('x', 'y')
return(out[,1])
}
以及一些基准测试数据:
set.seed(123456789)
sigma = rWishart(1, 2, diag(2))
means = rnorm(2)
X = rmvnorm(10000, means, sigma[,,1])
x1 = X[,1]
x2 = X[,2]
benchmark(mvnormDeriv(x1,x2,mu=means,sigma=sigma),
order="relative", replications=5)[,1:4]
公式可以在《矩阵手册》(2012)的第346个公式中找到。
我无法修改来自此处的多元正态分布的
rcpp
实现。以下是我尝试使用的一些代码。// [[Rcpp::export]]
arma::vec dmvnormDeriv_arma(arma::mat x, SEXP mu_sexp, arma::mat sigma, bool log = false) {
// create Rcpp vector and matrix from SEXP arguments
Rcpp::NumericVector mu_rcpp(mu_sexp);
// create views for arma objects(reuses memory and avoids extra copy)
arma::vec mu_vec(mu_rcpp.begin(), mu_rcpp.size(), false);
arma::rowvec mu(mu_rcpp.begin(), mu_rcpp.size(), false);
// return(mu_vec);
arma::vec distval = Mahalanobis(x, mu, sigma);
double logdet = sum(arma::log(arma::eig_sym(sigma)));
double log2pi = std::log(2.0 * M_PI);
arma::vec val = exp(-( (x.n_cols * log2pi + logdet + distval)/2));
// x.each_row() -= mu;
// arma::vec val2 = solve(sigma, x.row(1));
// arma::vec retval = -1 * val(1) * solve(sigma, x.row(1)-mu_vec);
return(val);
}
当然,这还不完整。你有什么想法可以在rcpp或者Armadillo中实现“* solve(sigma,(x-mu))”部分吗?我遇到了处理不同变量类型和运行每一行x的解决方案的问题。