r - Beta 分布的特征函数

Question

我正在尝试使用 R 计算许多不同 alpha 和 beta 的 beta 分布的特征函数；不幸的是，我遇到了数值问题。

首先我使用的是包 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 的其他值，我会得到错误:“检测到舍入错误”。

所以我的问题是:
对于 alpha 和 beta 的所有可能组合，是否有另一种方法可以获得 beta 分布特征函数的可靠结果？

我犯了什么明显的错误吗？

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")

如您所见 alpha=beta=1两种方法都提供相同的结果， cfX_Beta(t, alpha, beta)为 alpha=121.3618 疯狂和 beta=5041.483整合的结果似乎是合理的。对于 alpha=1和 beta=1/2集成失败。

Answer 1

它似乎适用于 RcppNumerical , 在不使用太小的公差的条件下 (1e-4 下面)。

// [[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")