- html - 出于某种原因,IE8 对我的 Sass 文件中继承的 html5 CSS 不友好?
- JMeter 在响应断言中使用 span 标签的问题
- html - 在 :hover and :active? 上具有不同效果的 CSS 动画
- html - 相对于居中的 html 内容固定的 CSS 重复背景?
我正在尝试让 Rcpp 版本的 pmvnorm 至少与 R 中的 mvtnorm::pmvnorm 一样快。
我找到了 https://github.com/zhanxw/libMvtnorm并使用相关源文件创建了一个 Rcpp 框架包。我添加了以下使用 Armadillo 的函数(因为我在编写的其他代码中使用了它)。
//[[Rcpp::export]]
arma::vec triangl(const arma::mat& X){
arma::mat LL = arma::trimatl(X, -1); // omit the main diagonal
return LL.elem(arma::find(LL != 0));
}
//[[Rcpp::export]]
double pmvnorm_cpp(arma::vec& bound, arma::vec& lowtrivec){
double error;
int n = bound.n_elem;
double* boundptr = bound.memptr();
double* lowtrivecptr = lowtrivec.memptr();
double result = pmvnorm_P(n, boundptr, lowtrivecptr, &error);
return result;
}
从构建包后的 R 中,这是一个可重现的示例:
set.seed(1)
covar <- rWishart(1, 10, diag(5))[,,1]
sds <- diag(covar) ^-.5
corrmat <- diag(sds) %*% covar %*% diag(sds)
triang <- triangl(corrmat)
bounds <- c(0.5, 0.9, 1, 4, -1)
rbenchmark::benchmark(pmvnorm_cpp(bounds, triang),
mvtnorm::pmvnorm(upper=bounds, corr = corrmat),
replications=1000)
这表明 pmvnorm_cpp 比 mvtnorm::pmvnorm 慢得多。结果是不同的。
> pmvnorm_cpp(bounds, triang)
[1] 0.04300643
> mvtnorm::pmvnorm(upper=bounds, corr = corrmat)
[1] 0.04895361
这让我感到困惑,因为我认为基本的 Fortran 代码是相同的。我的代码中有什么东西让一切变慢了吗?或者我应该尝试直接移植 mvtnorm::pmvnorm 代码吗?我几乎没有使用 Fortran 的经验。
感谢建议,请原谅我的无能。
编辑:为了与替代品进行快速比较,这个:
//[[Rcpp::export]]
NumericVector pmvnorm_cpp(NumericVector bound, NumericMatrix cormat){
Environment stats("package:mvtnorm");
Function f = stats["pmvnorm"];
NumericVector lower(bound.length(), R_NegInf);
NumericVector mean(bound.length());
NumericVector res = f(lower, bound, mean, cormat);
return res;
}
具有与 R 调用基本相同的性能(以下是 40 维 mvnormal):
> rbenchmark::benchmark(pmvnorm_cpp(bounds, corrmat),
+ mvtnorm::pmvnorm(upper=bounds, corr = corrmat),
+ replications=100)
test replications elapsed relative user.self sys.self
2 mvtnorm::pmvnorm(upper = bounds, corr = corrmat) 100 16.86 1.032 16.60 0.00
1 pmvnorm_cpp(bounds, corrmat) 100 16.34 1.000 16.26 0.01
所以在我看来,前面的代码中一定有什么事情发生了。无论是我如何处理 Armadillo 的事情,还是其他事情之间的联系。我认为与上一个实现相比应该有性能提升。
最佳答案
我不会为此尝试使用额外的库,而是尝试使用 mvtnorm
导出的 C API,c.f. https://github.com/cran/mvtnorm/blob/master/inst/NEWS#L44-L48 .在这样做的同时,我发现了导致结果不同的三个原因。其中之一还对性能差异负责:
mvtnorm
使用 R 的 RNG,而这已从您正在使用的库中删除,c.f. https://github.com/zhanxw/libMvtnorm/blob/master/libMvtnorm/randomF77.c .
您的 triangl
函数不正确。它以列优先顺序返回下三角矩阵。然而,底层的 Fortran 代码期望它按行优先顺序排列,c.f. https://github.com/cran/mvtnorm/blob/master/src/mvt.f#L36-L39和 https://github.com/zhanxw/libMvtnorm/blob/master/libMvtnorm/mvtnorm.cpp#L60
libMvtnorm
使用 1e-6
而不是 1e-3
作为相对精度,c.f. https://github.com/zhanxw/libMvtnorm/blob/master/libMvtnorm/mvtnorm.cpp#L65 .这也是造成性能差异的原因。
我们可以使用以下代码进行测试:
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::depends(mvtnorm)]]
#include <mvtnormAPI.h>
//[[Rcpp::export]]
arma::vec triangl(const arma::mat& X){
int n = X.n_cols;
arma::vec res(n * (n-1) / 2);
for (int i = 0; i < n; ++i) {
for (int j = 0; j < i; ++j) {
res(j + i * (i-1) / 2) = X(i, j);
}
}
return res;
}
// [[Rcpp::export]]
double pmvnorm_cpp(arma::vec& bound,
arma::vec& lowertrivec,
double abseps = 1e-3){
int n = bound.n_elem;
int nu = 0;
int maxpts = 25000; // default in mvtnorm: 25000
double releps = 0; // default in mvtnorm: 0
int rnd = 1; // Get/PutRNGstate
double* bound_ = bound.memptr();
double* correlationMatrix = lowertrivec.memptr();
double* lower = new double[n];
int* infin = new int[n];
double* delta = new double[n];
for (int i = 0; i < n; ++i) {
infin[i] = 0; // (-inf, bound]
lower[i] = 0.0;
delta[i] = 0.0;
}
// return values
double error;
double value;
int inform;
mvtnorm_C_mvtdst(&n, &nu, lower, bound_,
infin, correlationMatrix, delta,
&maxpts, &abseps, &releps,
&error, &value, &inform, &rnd);
delete[] (lower);
delete[] (infin);
delete[] (delta);
return value;
}
/*** R
set.seed(1)
covar <- rWishart(1, 10, diag(5))[,,1]
sds <- diag(covar) ^-.5
corrmat <- diag(sds) %*% covar %*% diag(sds)
triang <- triangl(corrmat)
bounds <- c(0.5, 0.9, 1, 4, -1)
set.seed(1)
system.time(cat(mvtnorm::pmvnorm(upper=bounds, corr = corrmat), "\n"))
set.seed(1)
system.time(cat(pmvnorm_cpp(bounds, triang, 1e-6), "\n"))
set.seed(1)
system.time(cat(pmvnorm_cpp(bounds, triang, 0.001), "\n"))
*/
结果:
> system.time(cat(mvtnorm::pmvnorm(upper=bounds, corr = corrmat), "\n"))
0.04896221
user system elapsed
0.000 0.003 0.003
> system.time(cat(pmvnorm_cpp(bounds, triang, 1e-6), "\n"))
0.04895756
user system elapsed
0.035 0.000 0.035
> system.time(cat(pmvnorm_cpp(bounds, triang, 0.001), "\n"))
0.04896221
user system elapsed
0.004 0.000 0.004
在相同的 RNG(和 RNG 状态)、正确的下三角相关矩阵和相同的相对精度下,结果相同且性能相当。精度越高,性能就越差。
所有这些都是针对使用 Rcpp::sourceCpp
的独立文件。为了在包中使用它,您需要将 LinkingTo: mvtnorm
添加到您的 DESCRIPTION
文件。
关于mvtnorm::pmvnorm 的 Rcpp 实现比原始 R 函数慢,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/51290014/
C语言sscanf()函数:从字符串中读取指定格式的数据 头文件: ?
最近,我有一个关于工作预评估的问题,即使查询了每个功能的工作原理,我也不知道如何解决。这是一个伪代码。 下面是一个名为foo()的函数,该函数将被传递一个值并返回一个值。如果将以下值传递给foo函数,
CStr 函数 返回表达式,该表达式已被转换为 String 子类型的 Variant。 CStr(expression) expression 参数是任意有效的表达式。 说明 通常,可以
CSng 函数 返回表达式,该表达式已被转换为 Single 子类型的 Variant。 CSng(expression) expression 参数是任意有效的表达式。 说明 通常,可
CreateObject 函数 创建并返回对 Automation 对象的引用。 CreateObject(servername.typename [, location]) 参数 serv
Cos 函数 返回某个角的余弦值。 Cos(number) number 参数可以是任何将某个角表示为弧度的有效数值表达式。 说明 Cos 函数取某个角并返回直角三角形两边的比值。此比值是
CLng 函数 返回表达式,此表达式已被转换为 Long 子类型的 Variant。 CLng(expression) expression 参数是任意有效的表达式。 说明 通常,您可以使
CInt 函数 返回表达式,此表达式已被转换为 Integer 子类型的 Variant。 CInt(expression) expression 参数是任意有效的表达式。 说明 通常,可
Chr 函数 返回与指定的 ANSI 字符代码相对应的字符。 Chr(charcode) charcode 参数是可以标识字符的数字。 说明 从 0 到 31 的数字表示标准的不可打印的
CDbl 函数 返回表达式,此表达式已被转换为 Double 子类型的 Variant。 CDbl(expression) expression 参数是任意有效的表达式。 说明 通常,您可
CDate 函数 返回表达式,此表达式已被转换为 Date 子类型的 Variant。 CDate(date) date 参数是任意有效的日期表达式。 说明 IsDate 函数用于判断 d
CCur 函数 返回表达式,此表达式已被转换为 Currency 子类型的 Variant。 CCur(expression) expression 参数是任意有效的表达式。 说明 通常,
CByte 函数 返回表达式,此表达式已被转换为 Byte 子类型的 Variant。 CByte(expression) expression 参数是任意有效的表达式。 说明 通常,可以
CBool 函数 返回表达式,此表达式已转换为 Boolean 子类型的 Variant。 CBool(expression) expression 是任意有效的表达式。 说明 如果 ex
Atn 函数 返回数值的反正切值。 Atn(number) number 参数可以是任意有效的数值表达式。 说明 Atn 函数计算直角三角形两个边的比值 (number) 并返回对应角的弧
Asc 函数 返回与字符串的第一个字母对应的 ANSI 字符代码。 Asc(string) string 参数是任意有效的字符串表达式。如果 string 参数未包含字符,则将发生运行时错误。
Array 函数 返回包含数组的 Variant。 Array(arglist) arglist 参数是赋给包含在 Variant 中的数组元素的值的列表(用逗号分隔)。如果没有指定此参数,则
Abs 函数 返回数字的绝对值。 Abs(number) number 参数可以是任意有效的数值表达式。如果 number 包含 Null,则返回 Null;如果是未初始化变量,则返回 0。
FormatPercent 函数 返回表达式,此表达式已被格式化为尾随有 % 符号的百分比(乘以 100 )。 FormatPercent(expression[,NumDigitsAfterD
FormatNumber 函数 返回表达式,此表达式已被格式化为数值。 FormatNumber( expression [,NumDigitsAfterDecimal [,Inc
我是一名优秀的程序员,十分优秀!