R/dChern.R

In ChernoffDist: Chernoff's Distribution

#' Density function of Chernoff's distribution
#'
#' Computes the density of Chernoff's distribution.
#'
#' @param x evaluation point of the density.
#' @return The function returns Chernoff's density evaluated at x.
#' @export
#' @examples
#' dChern(0)
#' @importFrom gsl airy_Ai_deriv airy_zero_Ai
dChern <- function(x){

  m1=20
  m2=100

  a <- c(1.0000000000000000000000000000000000000000000000000,
         0.14583333333333333333333333333333333333333333333333,
         -0.0014105902777777777777777777777777777777777777777778,
         -0.000020452697861552028218694885361552028218694885361552,
         1.7621771897199821461974239752017529795307573085351*10^-6,
         -5.7227597511835045485115756279777443798607819771841*10^-8,
         1.2738111577329144359096636472671763331236997246611*10^-9,
         -2.1640824986157345585013088038402781069395040086197*10^-11,
         2.8792575963578932547707953671069741538768134987350*10^-13,
         -2.9044965825655693838957326026543396338028657917270*10^-15,
         1.8213982724418504633090392507915213304919452526115*10^-17,
         4.3188090863188565193332245123091487021119622809634*10^-20,
         -3.4596348303488590277613451243743754391002537464995*10^-21,
         6.5412953607524627299042478809831523007851994829431*10^-23,
         -8.9227552497408589635214480311601977115813439133660*10^-25,
         1.0102911018761961214130852526736246207478858671070*10^-26,
         -9.9767598824744308308158851597119799263659899008997*10^-29,
         8.7968070582324563459058672533714015023352054588072*10^-31,
         -7.0120709194869959645933330651103032439601730987935*10^-33,
         5.0837836130061037314804057915251548236233281987597*10^-35,
         -3.3556860429372658821846021986263568219938365369527*10^-37,
         2.0068536445382019503549952854886809387341804072727*10^-39,
         -1.0717059663091867051094477160968394469729115856303*10^-41,
         4.9280004390027629026693452692753922980490376905806*10^-44,
         -1.7520728636583853900751602535859412088801567453800*10^-46,
         2.5147301934257098108652520310092365947204543091883*10^-49,
         3.0777765515961761642421078692169992079465879474309*10^-51,
         -4.1294815651393201362478292405036276929226139031707*10^-53,
         3.4347447301437443635331559664892799016709792483777*10^-55,
         -2.3612230361831206189177703747563180622623952194477*10^-57)

  b <- c(2/3, 4/189, -(64/81081), 976/49997493, -(77696/243952053345),
  572288/229236556407705, 18533023232/511729578836604433305,
  -(38925099776/20766274238875709829285),
  78619597150208/1789255032641387391955709835,
  -(179894925477619712/237567505371163930861834656295575),
  355939707144347451392/33483862154131262497617716288736895725,
  -(39005836704933300156915712/
    308463990815431083031137573116726428713954225),
  2289050102982895009959182336/
  1767060323806524493176742466881255509500679984825,
  -(200649578523736545985969979392/
    17542687184105536785979806459887088345436232597132175),
  4150584049080095253689757917642752/
  48582269248125992981340700593429828317979301804913004911625,
  -(82316715721780220361186959715104915456/
    162719863876102256856570365532382545656620198269910211176772324625),
  2496269073906379289258378706435696965976064/
  1395400775552373505526125181645270340244887721239451347831976357602870375, 40443356788234234398215577784896849409736704/
  5751763011877324021467054444486239515866394457455363392479567000207469259125, -(12400981672806616511218129306850440182537989914624/
            56378311862895846922417745137984433785556175656086728114156557759925453897236541375),
  2176427212822939497732275756764259696648526229405696/
  782306410974150617587789350479504850402170468932260389934546804909907969895916913145625,
  -(21561614139657230020572049244849602960430362030533246976/
      797222499419309852075034430608056068383121510467474191190315332840779470756715748947280724375),
  302827102989233974106114471447525975810838287995351689330688/
  1339244203480467350586528487628986627967035849852526411984858293733949481732299722715338387939360625,
  -(91775421041340022653388492161953065378621119513255608445829125413160453926836820815234058047522630993689775152134338585774820988788657566339842811063183943612385625),
  3604717664110903209581858970835165280561809951297369749960874000384/
309731916491210695594687805646165143740444845049895549634303721396184814320941371480033155514085647811463570625,
  -(670330408914808139855636443316181612197328948598502322269343700930789376/9058034348311231571679408604299235935704924155248341962370470451391273528114599422558499165462718158814590002422046875), 12481554650074942180082900504407664725598130549774321620122475633513070592/28430794314739805157391956315310763972644555715732930228398525551063595984630822795103411402181531931198755900576614859375,
  -(2178604749605281059085670479107218647441100172597386489603474027552121470083661824892788945832550571891557717372190617887147058341195226814720723678590396350063969327480262722082772114269771539265447416225817765625),
  118650261618618705922938600993055699827644959397645322131054710228993429756159459328/
  9315327591336862815891023166758399357976014558043175290733224414450712969357663236783314814967655424381057750345338837421114053922515625,
  -(416875163385443955693262949517093691729596823498254757123190231429866091828009095921664/
      6675191093438268092014849872334792085872395593140186845648574792607978187836199660185906738221925848707437339641027363569146502485038040015625), 12900593393454343659470432588607442580218696269127678503728505165469347629315460594292228096/
  44933776212674930334380945354056642398305849517717762531492550513413963261299009646409242716531978831467821737920209460399793121393353419025006015625)

  # Define the function p
  p <- function(t){

    sum_a <- 0
    sum_b <- 0
    p2 <- 0

    for (k in 1:m1){

      sum_a <- sum_a + a[k]*t^(3*(k-1))
      sum_b <- sum_b + b[k]*t^(3*(k-1/2))
      p2 <- p2 + exp(2^(1/3)*t*gsl::airy_zero_Ai(k))

    }

    return(
      (t<=1) * (-(pi/2)^(1/2)*sum_a + sum_b) +
      (t>1) * (-t^(-3/2) + 2*(2*pi)^(1/2)*exp(-t^3/6)*p2)
    )

  }

  h <- function(x){
    sapply(x, function(x_elem) {
      return(2*(2/pi)^(1/2)*stats::integrate(function(v){((2*x_elem + v^2)*v^2 + (1/2)*(2*x_elem + v^2)^2)*
          exp(-v^2*(2*x_elem + v^2)^2/2)},lower=0,upper=Inf)$value)
    })
  }

  g1 <- function(x){
    sapply(x, function(x_elem) {
      return((2*pi)^(-1/2)*stats::integrate(function(t){p(t)*exp(-t*(2*x_elem + t)^2/2)}, lower=0, upper=Inf)$value)
    })
  }

  g2 <- function(x){
    return(2*x - g1(x) + h(x))
  }

  g3 <- function(x){
    g3sum <- 0
    for (k in 1:m2){
      g3sum <- g3sum + exp(-2^(1/3)*gsl::airy_zero_Ai(k)*x)/gsl::airy_Ai_deriv(gsl::airy_zero_Ai(k))
    }
    return(4^(1/3)*exp((2/3)*x^3) * g3sum)
  }

  g <- function(x){return((x>-1)*g2(x) + (x<=-1)*g3(x))}

  return( g(x)*g(-x)/2 )

}

dChern(0)