Rで数値積分の記事

　https://arupaka-arupaka.hatenadiary.org/entry/2020/06/18/155816 を利用して、色々な積分をする。

マニュアル

https://www.gnu.org/software/gsl/doc/html/integration.html

コード例

ベータ分布の積分

$\int^{1} _ 0 Binom(n| size=m,prob=1-(1-x)^{5.5}) \cdot \frac{(x)^{\alpha-1}(1-x)^{\beta -1}}{Beta(\alpha,\beta)} dx$

　Rでの実装。解析解。

諸事情でαとβが二乗されているが、必要ない。

d_suzuki3_ver5<-function(m,n,alpha1,beta1){
    gamma1<-5.5
  
    ans<-mapply(function(m,n,alpha1,beta1){v<-0:m;1.0/base::beta(alpha1^2,beta1^2)*sum(choose(m,v)*(-1)^(v)*(base::beta(gamma1*((v)+n-m)+beta1^2,alpha1^2)))},m,n,alpha1,beta1)
  
    choose(n,m)*ans
  
}

Gauss-Legendre integration　（なめらかな関数では速い）

The fixed-order Gauss-Legendre integration routines are provided for fast integration of smooth functions with known polynomial order. The n-point Gauss-Legendre rule is exact for polynomials of order 2n - 1 or less. For example, these rules are useful when integrating basis functions to form mass matrices for the Galerkin method. Unlike other numerical integration routines within the library, these routines do not accept absolute or relative error bounds.

// [[Rcpp::depends(RcppGSL)]]
#include <RcppGSL.h>
#include <gsl/gsl_integration.h>
#include<gsl/gsl_sf.h>
#include <gsl/gsl_randist.h>
using namespace std;

struct my_f_params { double x; int size; double p1;double p2;};

double integrand(double p0, void *p) {
  struct my_f_params * params = (struct my_f_params *)p;
  double x = params->x;
  int size = params->size;
  double p1 = params->p1;
  double p2 = params->p2;
  int  m = size;
  double alpha = p1*p1;
  double beta = p2*p2;

  double res =gsl_ran_binomial_pdf(x,1-pow(1-p0,5.5),m)*gsl_ran_beta_pdf(p0, alpha, beta);
    //printf("res:%lf\n",res);
  
  return(res);

}

//double x; int size; double p1;double p2;
// [[Rcpp::export]]
Rcpp::NumericVector d_suzuki3_ver4(Rcpp::NumericVector x, Rcpp::IntegerVector size, Rcpp::NumericVector p1, Rcpp::NumericVector p2) {
  Rcpp::NumericVector resVec(x.size());

 gsl_integration_glfixed_table *table1 = gsl_integration_glfixed_table_alloc(100);
  
  gsl_function F;
  F.function = &integrand;
  struct my_f_params params = {x[0],(size)[0],p1[0],p2[0]};
  
  for (int i=0; i < x.size(); ++i) {
    //mu = xVec[i];
    params.x = x[i];
    params.size = size[i];
    params.p1 = p1[i];
    params.p2 = p2[i];
    F.params = &params;
  
    try{

      double result2=gsl_integration_glfixed(&F, 0,1,table1);
      resVec[i] = result2;
      
 
    }catch (char *str){

      resVec[i] = std::numeric_limits<double>::infinity();
      
    }
    
    gsl_integration_glfixed_table_free(table1);
  }

  return resVec;
}

QAG adaptive integration （積分区切りを自動調整）

The QAG algorithm is a simple adaptive integration procedure. The integration region is divided into subintervals, and on each iteration the subinterval with the largest estimated error is bisected. This reduces the overall error rapidly, as the subintervals become concentrated around local difficulties in the integrand. These subintervals are managed by the following struct,

// [[Rcpp::depends(RcppGSL)]]
#include <RcppGSL.h>
#include <gsl/gsl_integration.h>
#include<gsl/gsl_sf.h>
#include <gsl/gsl_randist.h>
using namespace std;

struct my_f_params { double x; int size; double p1;double p2;};

double integrand(double p0, void *p) {
  struct my_f_params * params = (struct my_f_params *)p;
  double x = params->x;
  int size = params->size;
  double p1 = params->p1;
  double p2 = params->p2;
  int  m = size;
  double alpha = p1*p1;
  double beta = p2*p2;
  double res =gsl_ran_binomial_pdf(x,1-pow(1-p0,5.5),m)*gsl_ran_beta_pdf(p0, alpha, beta);
  

  return(res);
    
}

//double x; int size; double p1;double p2;
// [[Rcpp::export]]
Rcpp::NumericVector d_suzuki3_ver4(Rcpp::NumericVector x, Rcpp::IntegerVector size, Rcpp::NumericVector p1, Rcpp::NumericVector p2) {
  Rcpp::NumericVector resVec(x.size());
 gsl_integration_workspace *work=gsl_integration_workspace_alloc(3000);
 // gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
    
 //size_t nevals;
  double result(0.0);
  double abserr(0.0);
  size_t nevals(300);
  //double x(x[0]);
  
  gsl_function F;
  F.function = &integrand;
  struct my_f_params params = {x[0],(size)[0],p1[0],p2[0]};
  
  for (int i=0; i < x.size(); ++i) {
    //mu = xVec[i];
    params.x = x[i];
    params.size = size[i];
    params.p1 = p1[i];
    params.p2 = p2[i];
    F.params = &params;
 
  
    try{
 
      int res2 = gsl_integration_qag(&F, 0,1, 1e-8, 1e-8,nevals,1,work, &result, &abserr);
      resVec[i] = result;

     
    }catch (char *str){
 
      resVec[i] = std::numeric_limits<double>::infinity();
      
    }
    
  }  
     gsl_integration_workspace_free (work);
    

  
 
  return resVec;
}

QAG adaptive integration （積分区切りを自動調整の特異点を配慮。両端の0と１が特異点）

ptsで特異点を指定している。

// [[Rcpp::depends(RcppGSL)]]
#include <RcppGSL.h>
#include <gsl/gsl_integration.h>
#include<gsl/gsl_sf.h>
#include <gsl/gsl_randist.h>
using namespace std;

struct my_f_params { double x; int size; double p1;double p2;};

double integrand(double p0, void *p) {
  struct my_f_params * params = (struct my_f_params *)p;
  double x = params->x;
  int size = params->size;
  double p1 = params->p1;
  double p2 = params->p2;
  int  m = size;
  double alpha = p1*p1;
  double beta = p2*p2;

  double res =gsl_ran_binomial_pdf(x,1-pow(1-p0,5.5),m)*gsl_ran_beta_pdf(p0, alpha, beta);
   
  return(res);
  
  
}

//double x; int size; double p1;double p2;
// [[Rcpp::export]]
Rcpp::NumericVector d_suzuki3_ver4(Rcpp::NumericVector x, Rcpp::IntegerVector size, Rcpp::NumericVector p1, Rcpp::NumericVector p2) {
  Rcpp::NumericVector resVec(x.size());
  //gsl_integration_cquad_workspace *work = gsl_integration_cquad_workspace_alloc(1e2);
 // gsl_integration_qaws_table *t = gsl_integration_qaws_table_alloc(alpha-1, double beta-1,0,0)
 //gsl_integration_romberg_workspace* work = gsl_integration_romberg_alloc(1e2);
 gsl_integration_workspace *work=gsl_integration_workspace_alloc(3000);
 //size_t nevals;
  double result(0.0);
  double abserr(0.0);
  size_t nevals(1000);
  double pts[2];
  pts[0]=0;
  pts[1]=1;
  //double x(x[0]);
  
  gsl_function F;
  F.function = &integrand;
  struct my_f_params params = {x[0],(size)[0],p1[0],p2[0]};
  
  for (int i=0; i < x.size(); ++i) {
    //mu = xVec[i];
    params.x = x[i];
    params.size = size[i];
    params.p1 = p1[i];
    params.p2 = p2[i];
    F.params = &params;
    //int res = gsl_integration_cquad(&F, 0.001,0.999, 1e-8, 1e-6, work, &result, &abserr, &nevals);
    //int res = gsl_integration_cquad(&F, 0.001,0.999, 1e-4, 1e-4, work, &result, &abserr, &nevals);
    //int res = gsl_integration_cquad(&F, 0,1, 1e-20, 1e-20, work, &result, &abserr, &nevals);
    //gsl_integration_qaws_table *t = gsl_integration_qaws_table_alloc(p1[i]*p1[i]-1,p2[i]*p2[i]-1,0,0);
    //gsl_integration_workspace *w = gsl_integration_workspace_alloc(1000);
    
    //int res = gsl_integration_qaws(&F, 0,1, t,1e-20, 1e-20,nevals,work, &result, &abserr);
    try{
 
      
      int res2 = gsl_integration_qagp(&F,pts,2, 1e-8, 1e-8,nevals,work, &result, &abserr);
      resVec[i] = result;
      

    }catch (char *str){
      //cout << str; 
      resVec[i] = std::numeric_limits<double>::infinity();
      
    }
    
  }
  gsl_integration_workspace_free (work);

  return resVec;
}

QAWS adaptive integration for singular functions (重み付き積分。重みが発散する)

The QAWS algorithm is designed for integrands with algebraic-logarithmic singularities at the end-points of an integration region. In order to work efficiently the algorithm requires a precomputed table of Chebyshev moments. 積分核が $(x-a)^-(-\alpha)*(x-b)^{-\beta}$ など $\alpha>0$ と $\beta>0$ など．これは、積分核が大きくなりすぎてうまくいかなかった。重み関数が端で発散する場合の積分。　重みと積分核でちいさくなるため。

// [[Rcpp::depends(RcppGSL)]]
#include <RcppGSL.h>
#include <gsl/gsl_integration.h>

#include <gsl/gsl_randist.h>

struct my_f_params { double x; int size; double p1;double p2;};

double integrand(double p0, void *p) {
  struct my_f_params * params = (struct my_f_params *)p;
  double x = params->x;
  int size = params->size;
  double p1 = params->p1;
  double p2 = params->p2;
  int  m = size;
  double alpha = p1*p1;
  double beta = p2*p2;
  double res =gsl_ran_binomial_pdf(x,p0,m)*gsl_ran_beta_pdf(1-pow(1-p0,1.0/5.5), alpha, beta)*1.0/5.5*pow(1-p0,1.0/5.5-1);
  return(res);
 
}


//double x; int size; double p1;double p2;
// [[Rcpp::export]]
Rcpp::NumericVector d_suzuki3_ver4(Rcpp::NumericVector x, Rcpp::IntegerVector size, Rcpp::NumericVector p1, Rcpp::NumericVector p2) {
  Rcpp::NumericVector resVec(x.size());
  gsl_integration_cquad_workspace *work = gsl_integration_cquad_workspace_alloc(1e2);
 
 //gsl_integration_romberg_workspace* work = gsl_integration_romberg_alloc(1e2);
 
  double result(0.0);
  double abserr(0.0);
  size_t nevals(1e2);
  //double x(x[0]);
  
  gsl_function F;
  F.function = &integrand;
  struct my_f_params params = {x[0],(size)[0],p1[0],p2[0]};
  
  for (int i=0; i < x.size(); ++i) {
    //mu = xVec[i];
    params.x = x[i];
    params.size = size[i];
    params.p1 = p1[i];
    params.p2 = p2[i];
    F.params = &params;
    //int res = gsl_integration_cquad(&F, 0.001,0.999, 1e-8, 1e-6, work, &result, &abserr, &nevals);
    //int res = gsl_integration_cquad(&F, 0.001,0.999, 1e-4, 1e-4, work, &result, &abserr, &nevals);
    int res = gsl_integration_cquad(&F, 0,1, 1e-20, 1e-20, work, &result, &abserr, &nevals);
    //int res=gsl_integration_romberg(&F, -0.0001,1.0001, 1e-20, 1e-20, work, &result, &abserr, &nevals);   
    resVec[i] = result;
  }
  //gsl_integration_romberg_workspacace_free(work);
  gsl_integration_cquad_workspace_free(work);
  return resVec;
}

これが正解。ベータ関数の逆数が大きくなりすぎて積分核が計算不能になるので、ベータ関数をくくりだしたバージョン

// [[Rcpp::depends(RcppGSL)]]
#include <RcppGSL.h>
#include <gsl/gsl_integration.h>
#include<gsl/gsl_sf.h>
#include <gsl/gsl_randist.h>
using namespace std;

struct my_f_params { double x; int size; double p1;double p2;};

double integrand(double p0, void *p) {
  struct my_f_params * params = (struct my_f_params *)p;
  double x = params->x;
  int size = params->size;
  double p1 = params->p1;
  double p2 = params->p2;
  int  m = size;
  double alpha = p1*p1;
  double beta = p2*p2;
  //double res =gsl_ran_binomial_pdf(x,1-p0,m)*gsl_ran_beta_pdf(1-pow(1-p0,1.0/5.5), alpha, beta)*1.0/5.5*pow(1-p0,1.0/5.5-1);
  double res =gsl_ran_binomial_pdf(x,1-pow(1-p0,5.5),m);
  //printf("res:%lf\n",res);
  if(res>=pow(10,15)){
    //res=std::numeric_limits<double>::infinity();
    res=1000;
  }
  
  
  //printf("res:%lf\n",res);
  
  return(res);
  
  
  

}

//double x; int size; double p1;double p2;
// [[Rcpp::export]]
Rcpp::NumericVector d_suzuki3_ver4(Rcpp::NumericVector x, Rcpp::IntegerVector size, Rcpp::NumericVector p1, Rcpp::NumericVector p2) {
  Rcpp::NumericVector resVec(x.size());

 gsl_integration_workspace *work=gsl_integration_workspace_alloc(3000);
 //size_t nevals;
  double result(0.0);
  double abserr(0.0);
  size_t nevals(300);
  //double x(x[0]);
  
  gsl_function F;
  F.function = &integrand;
  struct my_f_params params = {x[0],(size)[0],p1[0],p2[0]};
  
  for (int i=0; i < x.size(); ++i) {
    //mu = xVec[i];
    params.x = x[i];
    params.size = size[i];
    params.p1 = p1[i];
    params.p2 = p2[i];
    F.params = &params;
   
    gsl_integration_qaws_table *t = gsl_integration_qaws_table_alloc(p1[i]*p1[i]-1,p2[i]*p2[i]-1,0,0);
   
    //int res = gsl_integration_qaws(&F, 0,1, t,1e-20, 1e-20,nevals,work, &result, &abserr);
    try{
      int res2 = gsl_integration_qaws(&F, 0,1, t,1e-8, 1e-8,nevals,work, &result, &abserr);
      resVec[i] = result/gsl_sf_beta(p1[i]*p1[i],p2[i]*p2[i]);
      
    //int res=gsl_integration_romberg(&F, -0.0001,1.0001, 1e-20, 1e-20, work, &result, &abserr, &nevals);   
     
    }catch (char *str){
      //cout << str; 
      resVec[i] = std::numeric_limits<double>::infinity();
      
    }
    
    
    gsl_integration_qaws_table_free(t);
    //gsl_integration_cquad_workspace_free(work);
  }

  return resVec;
}

CQUAD doubly-adaptive integration （遅いけど二重の意味で発散に強い）

CQUAD is a new doubly-adaptive general-purpose quadrature routine which can handle most types of singularities, non-numerical function values such as Inf or NaN, as well as some divergent integrals. It generally requires more function evaluations than the integration routines in QUADPACK, yet fails less often for difficult integrands.

// [[Rcpp::depends(RcppGSL)]]
#include <RcppGSL.h>
#include <gsl/gsl_integration.h>

#include <gsl/gsl_randist.h>

struct my_f_params { double x; int size; double p1;double p2;};

double integrand(double p0, void *p) {
  struct my_f_params * params = (struct my_f_params *)p;
  double x = params->x;
  int size = params->size;
  double p1 = params->p1;
  double p2 = params->p2;
  int  m = size;
  double alpha = p1*p1;
  double beta = p2*p2;
  //double res =gsl_ran_binomial_pdf(x,p0,m)*gsl_ran_beta_pdf(1-pow(1-p0,1.0/5.5), alpha, beta)*1.0/5.5*pow(1-p0,1.0/5.5-1);
  double res =gsl_ran_binomial_pdf(x,1-pow(1-p0,5.5),m)*gsl_ran_beta_pdf(p0, alpha, beta);
  return(res);
  


}

//double x; int size; double p1;double p2;
// [[Rcpp::export]]
Rcpp::NumericVector d_suzuki3_ver4(Rcpp::NumericVector x, Rcpp::IntegerVector size, Rcpp::NumericVector p1, Rcpp::NumericVector p2) {
  Rcpp::NumericVector resVec(x.size());
  gsl_integration_cquad_workspace *work = gsl_integration_cquad_workspace_alloc(1e2);
 
 //gsl_integration_romberg_workspace* work = gsl_integration_romberg_alloc(1e2);
 
  double result(0.0);
  double abserr(0.0);
  size_t nevals(1e2);
  //double x(x[0]);
  
  gsl_function F;
  F.function = &integrand;
  struct my_f_params params = {x[0],(size)[0],p1[0],p2[0]};
  
  for (int i=0; i < x.size(); ++i) {
    //mu = xVec[i];
    params.x = x[i];
    params.size = size[i];
    params.p1 = p1[i];
    params.p2 = p2[i];
    F.params = &params;
    //int res = gsl_integration_cquad(&F, 0.001,0.999, 1e-8, 1e-6, work, &result, &abserr, &nevals);
    //int res = gsl_integration_cquad(&F, 0.001,0.999, 1e-4, 1e-4, work, &result, &abserr, &nevals);
    int res = gsl_integration_cquad(&F, 0,1, 1e-20, 1e-20, work, &result, &abserr, &nevals);
    //int res=gsl_integration_romberg(&F, -0.0001,1.0001, 1e-20, 1e-20, work, &result, &abserr, &nevals);   
    resVec[i] = result;
  }
  //gsl_integration_romberg_workspacace_free(work);
  gsl_integration_cquad_workspace_free(work);
  return resVec;
}

arupaka-_-arupakaの日記

Ｒで数値積分の続き。色々な積分手法を使う（コード例）

Rで数値積分の記事

マニュアル

コード例

Rでの実装。解析解。

Gauss-Legendre integration　（なめらかな関数では速い）

QAG adaptive integration （積分区切りを自動調整）

QAG adaptive integration （積分区切りを自動調整の特異点を配慮。両端の0と１が特異点）

QAWS adaptive integration for singular functions (重み付き積分。重みが発散する)

これが正解。ベータ関数の逆数が大きくなりすぎて積分核が計算不能になるので、ベータ関数をくくりだしたバージョン

CQUAD doubly-adaptive integration （遅いけど二重の意味で発散に強い）

Rで数値積分の記事

マニュアル

コード例

Rでの実装。解析解。

Gauss-Legendre integration （なめらかな関数では速い）

QAG adaptive integration （積分区切りを自動調整）

QAG adaptive integration （積分区切りを自動調整の特異点を配慮。両端の0と１が特異点）

QAWS adaptive integration for singular functions (重み付き積分。重みが発散する)

これが正解。ベータ関数の逆数が大きくなりすぎて積分核が計算不能になるので、ベータ関数をくくりだしたバージョン

CQUAD doubly-adaptive integration （遅いけど二重の意味で発散に強い）

　Rでの実装。解析解。

Gauss-Legendre integration　（なめらかな関数では速い）