如何在 Rcpp 中以数值方式计算积分 [英] How to calculate integral, numerically, in Rcpp
本文介绍了如何在 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 ac++ 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:
- 我真的不知道如何选择步骤,或者我需要多少个子区间来近似积分.我在本科的时候学过数值分析,我想我可能需要查一下我的关于误差项表达的书,来决定步长.
- 我将我的结果与 R 的结果进行了比较.R 中的integr() 函数可以处理区间[0,1] 内的积分.这对我有帮助,因为我的函数在 0 或 1 处未定义,它具有无限值.在我的 C++ 代码中,我只能从 [1e-4, 1-1e-4] 设置间隔.我尝试了不同的值,如 1e-7、1e-10,但是,1e-4 是最接近 R 的结果....我该怎么办?
这篇关于如何在 Rcpp 中以数值方式计算积分的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!
查看全文