使用逆 CDF 在 R 中生成随机变量答案

【问题标题】：Use inverse CDF to generate random variable in R使用逆 CDF 在 R 中生成随机变量
【发布时间】：2014-06-09 02:57:46
【问题描述】：

首先，我不知道教授是否给出了错误的问题。无论如何，我尝试生成F(x)~U(0,1)，其中CDF F(x)=1-(1+x)exp(-x)（对于此CDF，您无法手动计算x=g(F(x))）。然后计算F(x)的根，就可以实现问题想要的了。

因为从0 到INF 的根范围，uniroot() 是没有问题的。因此，我用牛顿法来写一个。

那么，我的代码是这样的：

f=function(x) {
      ifelse(x>=0,x*exp(-x),0)
  }

in.C=function(n) {
         a=runif(n)
         G=NULL
         for(i in 1:n) {
             del=1
             x=2
             while(abs(del)>1e-12){
               del=(1-(1+x)*exp(-x)-a[i])/f(x)
               x=x-del

             }
             G[i]=x
         }
         G
     }
system.time(tt<-in.C(100000))

但是，如果F(x) 太小，而在牛顿法中一步，结果可能小于零，则会发生错误。此外，我这样修改了我的代码：

f=function(x) {
      ifelse(x>=0,x*exp(-x),0)
  }

in.C=function(n) {
         a=runif(n)
         G=NULL
         for(i in 1:n) {
             del=1
             x=2
             while(abs(del)>1e-12){
               if(x>=0){    del=(1-(1+x)*exp(-x)-a[i])/f(x)
                   x=x-del
                   }
                   else break
             }
             if(x>=0) G[i]=x
         }
         G[!is.na(G)]
     }
system.time(tt<-in.C(100000))
hist(tt, breaks=70, right=F, freq=F)
curve(f(x),from=0,to=20,add=T)

很明显，代码是错误的，因为我拒绝了接近零的结果。

所以，我的问题是我的代码是否可以修改以计算正确，如果不能，是否有其他方法可以做到。感谢您提供任何帮助。

【问题讨论】：

标签： r newtons-method cdf

【解决方案1】：

您可以为此使用uniroot(...)。

[注意：如果本练习的目的是实现您自己的牛顿拉夫森技术版本，请告诉我，我将删除答案。]

如果我理解正确，您希望从概率密度函数f 和累积密度F 的分布中生成随机样本

f = x*exp(-x)
F = 1 - (1+x)*exp(-x)

正如您所暗示的，这可以通过从 U[0,1] 生成随机样本并根据 F 的逆 CDF 进行转换来完成。该过程与 here 和 here 发布的过程非常相似，只是您已经有了 CDF 的表达式。

f <- function(x) x*exp(-x)
F <- function(x) 1-(1+x)*exp(-x)

F.inv <- function(y){uniroot(function(x){F(x)-y},interval=c(0,100))$root}
F.inv <- Vectorize(F.inv)

x <- seq(0,10,length.out=1000)
y <- seq(0,1,length.out=1000)

par(mfrow=c(1,3))
plot(x,f(x),type="l",main="f(x)")
plot(x,F(x),type="l",main="CDF of f(x)")
plot(y,F.inv(y),type="l",main="Inverse CDF of f(x)")

然后，生成X ~ U[0,1] 和Z = F.inv(X)。

set.seed(1)
X <- runif(1000,0,1)   # random sample from U[0,1]
Z <- F.inv(X)

par(mfrow=c(1,1))
hist(Z, freq=FALSE, breaks=c(seq(0,10,length=30),Inf), xlim=c(0,10))
lines(x,f(x),type="l",main="Density function", col="red",lty=2)

【讨论】：