ベータ分布の特性関数

Question

私は、多くの異なるアルファとベータのために、ベータ分布の特性関数をRで計算しようとしています。残念ながら私は数値的な問題にぶつかっています。

まず以下の例を参照してください（私はパッケージCharFun、ほとんどの場合のために正常に動作しているようだ機能cfX_Beta(t, alpha, beta)を使用していたが、alpha=121.3618とbeta=5041.483ため例えば、それは完全に失敗し、Re(cfX_Beta(t, alpha, beta))とIm(cfX_Beta(t, alpha, beta))は常に[間隔にする必要があります-1,1]）。

次に、統合によって特性関数を求めることにしました。この方法ではalpha=121.3618とbeta=5041.483の信憑性のある結果が得られますが、他の組み合わせの場合は統合が失敗します。 （エラー：「積分はおそらく発散している」）。積分のためにrel.tolを増やしても助けにならなかった。逆に、他のアルファ値とベータ値については、「丸め誤差が検出されました」というエラーが表示されます。

質問： アルファとベータのすべての可能な組み合わせについて、ベータ版の特性関数の信頼できる結果を得る別の方法はありますか？

私は明らかな間違いをしていますか？

alpha=beta=1両方の方法は、同じ結果を提供するためにあなたが見ることができるように

library(CharFun) 
abc<-function(x,t,a,b) { 
    return(dbeta(x,a,b)*cos(t*x)) 
} 
dfg<-function(x,t,a,b) { 
    return(dbeta(x,a,b)*sin(t*x)) 
} 
hij<-function(t,a,b) { 
    intRe=rep(0,length(t)) 
    intIm=rep(0,length(t)) 
    i<-complex(1,0,1) 
    for (j in 1:length(t)) { 
    intRe[j]<-integrate(abc,lower=0,upper=1,t[j],a,b)$value 
    intIm[j]<-integrate(dfg,lower=0,upper=1,t[j],a,b)$value 
    } 
    return(intRe+intIm*i) 
} 

alpha<-1 
beta<-1 

t <- seq(-100, 100, length.out = 501) 
par(mfrow=c(3,2)) 
alpha<-1 
beta<-1 
plotGraf(function(t) 
    hij(t, alpha, beta), t, title = "CF alpha=1 
beta=1") 
plotGraf(function(t) 
    cfX_Beta(t, alpha, beta), t, title = "CF Charfun alpha=1 
beta=1") 

alpha<-121.3618 
beta<-5041.483 
plotGraf(function(t) 
    hij(t, alpha, beta), t, title = "CF alpha=121.3618 beta=5041.483") 
plotGraf(function(t) 
    cfX_Beta(t, alpha, beta), t, title = "CF Charfun alpha=121.3618 beta=5041.483") 

alpha<-1 
beta<-1/2 
plotGraf(function(t) 
    hij(t, alpha, beta), t, title = "CF alpha=1 
beta=1/2") 
plotGraf(function(t) 
    cfX_Beta(t, alpha, beta), t, title = "CF Charfun alpha=1 
beta=1/2") 

、 cfX_Beta(t, alpha, beta)は alpha=121.3618と beta=5041.483統合の結果はもっともらしく思えるのために野生行きます。 alpha=1と beta=1/2の場合、統合は失敗します。

Answer 1

あまりにも小さな公差（1e-4）を使用する条件では、RcppNumericalと動作するようです。

// [[Rcpp::depends(RcppEigen)]] 
// [[Rcpp::depends(RcppNumerical)]] 
#include <RcppNumerical.h> 
using namespace Numer; 

class BetaCDF_Re: public Func 
{ 
private: 
    double a; 
    double b; 
    double t; 
public: 
    BetaCDF_Re(double a_, double b_, double t_) : a(a_), b(b_), t(t_){} 

    double operator()(const double& x) const 
    { 
    return R::dbeta(x, a, b, 0) * cos(t*x); 
    } 
}; 

class BetaCDF_Im: public Func 
{ 
private: 
    double a; 
    double b; 
    double t; 
public: 
    BetaCDF_Im(double a_, double b_, double t_) : a(a_), b(b_), t(t_) {} 

    double operator()(const double& x) const 
    { 
    return R::dbeta(x, a, b, 0) * sin(t*x); 
    } 
}; 

// [[Rcpp::export]] 
Rcpp::List integrate_test(double a, double b, double t) 
{ 
    BetaCDF_Re f1(a, b, t); 
    double err_est1; 
    int err_code1; 
    const double res1 = integrate(f1, 0, 1, err_est1, err_code1, 
           100, 1e-4, 1e-4, 
           Integrator<double>::GaussKronrod201); 
    BetaCDF_Im f2(a, b, t); 
    double err_est2; 
    int err_code2; 
    const double res2 = integrate(f2, 0, 1, err_est2, err_code2, 
           100, 1e-4, 1e-4, 
           Integrator<double>::GaussKronrod201); 
    return Rcpp::List::create(
    Rcpp::Named("realPart") = 
     Rcpp::List::create(
     Rcpp::Named("value") = res1, 
     Rcpp::Named("error_estimate") = err_est1, 
     Rcpp::Named("error_code") = err_code1 
    ), 
    Rcpp::Named("imPart") = 
     Rcpp::List::create(
     Rcpp::Named("value") = res2, 
     Rcpp::Named("error_estimate") = err_est2, 
     Rcpp::Named("error_code") = err_code2 
    ) 
); 
} 

> integrate_test(1, 0.5, 1) 
$realPart 
$realPart$value 
[1] 0.7497983 

$realPart$error_estimate 
[1] 7.110548e-07 

$realPart$error_code 
[1] 0 


$imPart 
$imPart$value 
[1] 0.5934922 

$imPart$error_estimate 
[1] 5.54721e-07 

$imPart$error_code 
[1] 0 

プロット：

t <- seq(-100, 100, length.out = 501) 
x <- lapply(t, function(t) integrate_test(1,0.5,t)) 
realparts <- unlist(purrr::transpose(purrr::transpose(x)$realPart)$value) 
imparts <- unlist(purrr::transpose(purrr::transpose(x)$imPart)$value) 
plot(t, realparts, type="l", col="blue", ylim=c(-1,1)) 
lines(t, imparts, type="l", col="red")