r - 多元正态/高斯与rcpp和RcppArmadillo的一阶导数

我正在尝试基于发布在here和here上的多元正态分布的rcpp实现，在R中实现多元正态分布的一阶导数。

这是一个快速的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中找到。

我无法从here修改多元正态的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);
}

当然，这还不完整。关于如何在* solve(sigma,(x-mu))或rcpp中实现Armadillo部分的任何想法？我在处理不同的变量类型并为x的每一行运行求解时遇到问题。

最佳答案

这是基于RcppArmadillo的解决方案。它比R实施快100倍以上。首先，c ++实现依赖于this rcpp gallery example。

// [[Rcpp::export]]
arma::mat dmvnormderiv_arma(arma::mat x, arma::rowvec mean, arma::mat sigma, bool log = false) {
    // get result for mv normal
    arma::vec distval = Mahalanobis(x,  mean, sigma);
    double logdet = sum(arma::log(arma::eig_sym(sigma)));
    double log2pi = std::log(2.0 * M_PI);
    arma::vec mvnorm = exp(-( (x.n_cols * log2pi + logdet + distval)/2));

    // create output matrix with one column for each derivative
    int n = x.n_rows;
    arma::mat deriv;
    deriv.copy_size(x);
    for (int i=0; i < n; i++) {
        deriv.row(i) = -1 * mvnorm(i) * trans(solve(sigma, trans(x.row(i) - mean)));
    }

    return(deriv);
}

还有两个R实现。一个是纯R，另一个是基于dmvnorm包中的mvtnorm。

library('RcppArmadillo')
library('mvtnorm')
library('rbenchmark')
sourceCpp('mvnorm.cpp')

mvnormDeriv = function(X, mu=rep(0,ncol(X)), sigma=diag(ncol(X))) {
    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(X,1,fn))
    return(out)
}
dmvnormDeriv = function(X, mean, sigma) {
    if (is.vector(X)) X <- matrix(X, ncol = length(X))
    if (missing(mean)) mean <- rep(0, length = ncol(X))
    if (missing(sigma)) sigma <- diag(ncol(X))
    n = nrow(X)
    mvnorm = dmvnorm(X, mean = mean, sigma = sigma)
    deriv = array(NA,c(n,ncol(X)))
    for (i in 1:n)
        deriv[i,] = -mvnorm[i] * solve(sigma,(X[i,]-mean))
    return(deriv)
}

最后是一些基准：

set.seed(123456789)
sigma = rWishart(1, 2, diag(2))[,,1]
means = rnorm(2)
X     = rmvnorm(10000, means, sigma)

benchmark(dmvnormderiv_arma(X,means,sigma),
        mvnormDeriv(X,mu=means,sigma=sigma),
        dmvnormDeriv(X,mean=means,sigma=sigma),
        order="relative", replications=5)[,1:4]

                                          test replications elapsed
1           dmvnormderiv_arma(X, means, sigma)            5   0.016
3 dmvnormDeriv(X, mean = means, sigma = sigma)            5   2.118
2    mvnormDeriv(X, mu = means, sigma = sigma)            5   5.939
  relative
1    1.000
3  132.375
2  371.187