r - R 和 Rcpp 中的多个多元正态密度值

Question

我有一个关于快速实施的问题。想象一下，您有一个矩阵 Ys，其中每一行都指代源自多元正态分布的观察值向量，例如，

    Ys = matrix(c(1.0,1.0,1.0,0.0,0.5,0.6,0.1,0.1,0.3), nrow = 3, ncol = 3) 

此外，还有一个矩阵 Sigs，其中每一行表示 Ys 中每个结果向量的方差协方差矩阵的对角元素，例如，

    Sigs = matrix(c(1.0,0.5,0.1,0.2,0.3,0.4,0.3,0.7,0.8), nrow = 3, ncol = 3) 

我想要做的是计算 Ys 中每一行的密度值，给定 Sigs 中相应行中的对角线元素。

可以在 R 中使用 for 循环，例如

    colSigs = ncol(Sigs)
    res = rep(0,3)
    means = rep(0,colSigs)
    for (i in 1:nrow(Ys) ) {
            sigma = diag(Sigs[i,],colSigs)
            res[i] = mvtnorm::dmvnorm(Ys[i,],means,sigma)
    }

但是，在我的例子中，Ys 和 Sigs 包含大约 100,000 行。所以我写了一个更快的 Rcpp 函数。不过，我想知道是否有一个花哨的技巧（一种更有效的方法），这样我就不必循环了？欢迎任何想法。

----

编辑：我被要求添加 Rcpp 函数。干得好：

此函数计算出现在多元正态密度中的二次形式：

double dmvnorm_distance( arma::rowvec y, arma::mat Sigma )
    {
    int n = Sigma.n_rows;
    double res=0;
    double fac=1;
    for (int ii=0; ii<n; ii++){
        for (int jj=ii; jj<n; jj++){
            if (ii==jj){ fac = 1; } else { fac = 2;}
            res += fac *y(0,ii) * Sigma(ii,jj) * y(0,jj);
        }
    }
    return res;
}

此函数计算密度值：

双 dmvnorm_rcpp( arma::rowvec y, arma::mat Sigma ) {

    int p = Sigma.n_rows;

    // inverse Sigma
    arma::mat Sigma1 = arma::inv(Sigma);
    // determinant Sigma
    double det_Sigma = arma::det(Sigma);
    // distance
    double dist = dmvnorm_distance( y, Sigma1);

    double pi1 = 3.14159265358979;
    double l1 = - p * std::log(2*pi1) - dist - std::log( det_Sigma );
    double ll = 0.5 * l1;

    return ll;
}

这个函数包含for循环并从R调用：

Rcpp::NumericVector mvnorm_loop( arma::mat Ys, arma::mat SIGs )
{
    int n = Ys.n_rows;
    Rcpp::NumericVector out(n);

    for (int ii=0; ii<n; ii++){

         // get yi and diagonal entries
         arma::rowvec yi = Ys.row(ii);
         arma::rowvec si = SIGs.row(ii);; 

         // make Sigma
         arma::mat Sigma = arma::diagmat(si);

         // compute likelihood value       
         out[ii] = dmvnorm_rcpp( yi, Sigma );
    }

    return out;
}

所以基本上问题是是否有另一种方法来实现 Rcpp 中的插入以使整个事情变得更快。

----

最好的，斯特凡

PS：我也在 R 中使用了 apply ，它比 Rcpp 循环函数慢。

标签： r