如何在 Rcpp 中以数值方式计算积分 [英] How to calculate integral, numerically, in Rcpp

问题描述

我已经搜索了一个小时来寻找进行数值积分的方法.我是 Rcpp 的新手，现在正在重写我的旧程序.我在 R 中所做的是:

I've searched for an hour for the methods doing numerical integration. I'm new to Rcpp and rewriting my old programs now. What I have done in R was:

   x=smpl.x(n,theta.true)   
   joint=function(theta){# the joint dist for
                         #all random variable
     d=c()
     for(i in 1:n){
       d[i]=den(x[i],theta)
     }
     return(prod(d)*dbeta(theta,a,b))   }  
 joint.vec=Vectorize(joint)##vectorize the function, as required when
                           ##using integrate()   
margin=integrate(joint.vec,0,1)$value # the
                                ##normalizeing constant at the donominator  
 area=integrate(joint.vec,0,theta.true)$value # the values at the
                                              ## numeritor

• R 中的 integrate() 函数会很慢，而且由于我正在对大小为 n 的样本的后验分布进行积分，因此积分值会很大且误差很大.
• 我试图在 Rcpp 的帮助下重写我的代码，但我不知道如何处理集成.我应该包含一个 c++ h 文件吗?或者有什么建议?
• The integrate() function in R will be slow, and since I am doing the integration for a posterior distribution of a sample of size n, the value of the integration will be huge with large error.
• I am trying to rewrite my code with the help of Rcpp, but I　don't know how to deal with the integrate. Should I include a c++ h file? Or any suggestions?

推荐答案

在阅读了您的 @utobi 建议后，我觉得自己编程可能更容易.我只是用辛普森公式来近似积分:

after reading your @utobi advice, I felt programming by my own maybe easier. I simply use Simpson formula to approximate the integral:

// [[Rcpp::export]]
double den_cpp (double x, double theta){
return(2*x/theta*(x<=theta)+2*(1-x)/(1-theta)*(theta<x));
}
// [[Rcpp::export]]
double joint_cpp ( double theta,int n,NumericVector x, double a, double b){
double val = 1.0;
NumericVector d(n);
for (int i = 0; i < n; i++){
double tmp = den_cpp(x[i],theta);
val = val*tmp;
}
val=val*R::dbeta(theta,a,b,0);
return(val);
}
// [[Rcpp::export]]
List Cov_rate_raw ( double theta_true,  int n, double a, double b,NumericVector x){
//This function is used to test, not used in the fanal one
int steps = 1000;
double s = 0;
double start = 1.0e-4;
std::cout<<start<<" ";
double end = 1-start;
std::cout<<end<<" ";
double h = (end-start)/steps;
std::cout<<"1st h ="<<h<<" ";
double area = 0;
double margin = 0;
for (int i = 0; i < steps ; i++){
double at_x = start+h*i;
double f_val = (joint_cpp(at_x,n,x,a,b)+4*joint_cpp(at_x+h/2,n,x,a,b)+joint_cpp(at_x+h,n,x,a,b))/6;
s = s + f_val;
}
margin = h*s;
s=0;
h=(theta_true-start)/steps;
std::cout<<"2nd h ="<<h<<" ";
for (int i = 0; i < steps ; i++){
double at_x = start+h*i;
double f_val = (joint_cpp(at_x,n,x,a,b)+4*joint_cpp(at_x+h/2,n,x,a,b)+joint_cpp(at_x+h,n,x,a,b))/6;
s = s + f_val;
}
area = h * s;
double r = area/margin;
int cover = (r>=0.025)&&(r<=0.975);
List ret;
ret["s"] = s;
ret["margin"] = margin;
ret["area"] = area;
ret["ratio"] = r;
ret["if_cover"] = cover;
return(ret);
}

我不太擅长 C++，所以两个 for 循环有点傻.

I'm not that good at c++, so the two for loops like kind of silly.

大体上可行，但仍有几个潜在问题:

It generally works, but there are still several potential problems:

1. 我真的不知道如何选择步骤，或者我需要多少个子区间来近似积分.我在本科的时候学过数值分析，我想我可能需要查一下我的关于误差项表达的书，来决定步长.
2. 我将我的结果与 R 的结果进行了比较.R 中的integr() 函数可以处理区间[0,1] 内的积分.这对我有帮助，因为我的函数在 0 或 1 处未定义，它具有无限值.在我的 C++ 代码中，我只能从 [1e-4, 1-1e-4] 设置间隔.我尝试了不同的值，如 1e-7、1e-10，但是，1e-4 是最接近 R 的结果....我该怎么办?

这篇关于如何在 Rcpp 中以数值方式计算积分的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！