первая производная от многомерного нормального/гауссова с rcpp и RcppArmadillo

Question

Я пытаюсь реализовать первую производную от многомерного нормального распределения в R на основе rcpp реализации многомерного нормального распределения, размещенного here и here.

Вот быстрое внедрение 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 многомерного нормали с here. Вот какой код, который я использовал, чтобы попробовать

// [[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.

Answer 1

Предложение представляет собой решение основано на RcppArmadillo. Это более чем в 100 раз быстрее, чем реализация R. Во-первых, реализация 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