0

首先，我不知道教授是否给出了错误的问题。无论如何，我尝试生成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)

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

所以，我的问题是我的代码是否可以修改以计算正确，如果不能，是否有其他方法可以做到。任何帮助表示赞赏。

4

1 回答 1

2

你可以用uniroot(...)这个。

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

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

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

正如您所暗示的，这可以通过生成随机样本U[0,1]并根据 的逆 CDF 进行转换来完成F。该过程与此处和此处发布的过程非常相似，只是您已经有了 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)

于 2014-04-24T03:10:28.747 回答