r - 多次复制 uniroot 函数

Question

我需要执行概率积分变换的蒙特卡罗方法才能找到 Pr[2.3<a<3.8]。我有一个分布，其中“xvalues”是从 0 到 5 的 200 个值的数据集：

A<-function(a,beta){
  A<-atan(a/beta)+atan((5-a)/beta)
  return(A)
}

cond_posterior<-function(a){
  xvalues<-flashes
  posterior <- (0.5)^3*exp(-8*0.5)*(0.5/A(a,0.5))*prod(1/((0.5)^2+(xvalues-a)^2))
  return(posterior)
}

这就是它的样子：

alphas_c<-seq(0,5,length.out = 1000)
posteriors<-NULL
for (i in alphas_c){
  posterior<-posterior_c(i)
  posteriors<-c(posteriors,posterior)
}

plot(posteriors~alphas_c,type="l",xlim=c(2,3.2))

我是如何做到这一点的，是使用积分来找到累积分布函数（cdf）：

vec.cond_posterior <- Vectorize(FUN=cond_posterior, vectorize.args = c("a"))


cdf_cond<-function(a){
  prob<-integrate(vec.cond_posterior,lower = 0,upper=a)
  return(prob)
}

然后我必须遵循以下逻辑：

我必须从标准均匀分布生成一个随机值 (u)，并通过求解从均匀分布生成的值处的 cdf 的倒数来求解 a 的值。然而，确实如此。找到逆 cdf 有点困难。因此，我们求解 F(a) - u = 0。为了做到这一点，我编写了以下函数并使用 uniroot 来求解表达式 F(a) - u = 0。

cdf_solve<-function(a){
  u<-runif(1,0,1)
  alpha<-vec.cond_posterior(a)-u
  return(alpha)
}

uniroot(cdf_solve,interval = c(2.8,3.2))

我认为这可以为 a 找到一个值。但是我需要重复很多次，比如 5000 次。但是，我被困在如何做到这一点上。我尝试过使用 for 循环，但我不知道如何重写代码以便可以多次执行 uniroot 函数。请帮忙！

标签： rrepeatintegral