2014-04-24 5 views
0

먼저, 교수가 잘못된 질문을 던졌습니다. 어쨌든 F(x)~U(0,1)을 생성하려고 시도했습니다. 여기서 CDF F(x)=1-(1+x)exp(-x) (이 CDF의 경우 수동으로 x=g(F(x))을 계산할 수 없습니다). 그리고 나서 질문의 목적을 달성하기 위해 F(x)의 근원을 계산하십시오.역 (逆) CDF를 사용하여 무작위 변수를 생성하십시오.

루트 범위가 0에서 INF까지이므로 uniroot()이 문제가됩니다. 따라서 저는 Newton Method를 사용하여 하나를 씁니다.

그럼, 내 코드는 다음과 같다 : 그러나

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)가 너무 작고, 뉴턴 방법에서 한 단계, 결과가 0 미만, 다음 오류가 일어날 수 있습니다합니다. 또한 다음과 같이 코드를 수정했습니다.

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) 

0에 가까운 결과를 거부했기 때문에 분명히 코드가 잘못되었습니다.

그래서 내 판단은 내 코드가 올바르게 계산되도록 수정 될 수 있는지 아닌지, 그렇지 않은지 여부입니다. 어떤 assitance도 감사하겠습니다.

답변

1

uniroot(...)을 사용하십시오.

[참고 :.이 연습의 포인트는 뉴턴 랩슨 기술의 자신의 버전을 구현하는 경우에, 저에게 알려 주시면 답변을 삭제됩니다]

것은 내가 제대로 당신이 이해하고있어 경우 확률 밀도 함수 f하고 알 수 있듯이

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

, 이는 F을 그 역에있어서 CDF를 U[0,1]에서 임의의 샘플을 생성하고, 변환에 의해 수행 될 수 누적 밀도 F와 분포로부터 임의의 샘플을 생성 할 . 이 절차는 CDF에 대한 표현식을 이미 가지고 있다는 점을 제외하면 herehere으로 게시 된 절차와 매우 유사합니다.

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)