比オブ制服分布

Question

私は私は次のようにアルゴリズムを作成する必要た運動有する：制服の

比を事実に基づいているが、その密度fの確率変数X（xについて）我々は、V）均一セットに分布

Af = {(u,v):0 < v ≤ f(u/v)} 

ランダム点は最小 - から除去によってAfの中に均一にサンプリングすることができる、（対についてX = U/Vを計算することによって所望の密度からUをXを生成することができますすなわち、Afを含む最小の可能な長方形。 それは（U-、U +）×（0、Vが+）で与えられる

v+ = max f(x), x 

u− = minx f(x), x 

u+ = maxx f(x) 

そして比オブ制服方法は、以下の簡単な手順で構成されてここで

1. 乱数を生成しUを（u-、u +）で一様に含む。

2. （0、v +）に乱数Vを一様に生成します。

3. X←U/Vを設定します。

4. F（X）≤V 2受け入れ、返すX.エルス

5. が再びしよう。

これまでの私のコード：

x <- cnorm(1, mean = 0, sd=1) 

myrnorm <- function(pdf){ 
    ## call rou() n times 
    pdf <- function(x) {exp(-x^2/2)} 
    } 
rou <- function(u, v) { 
    uplus <- 1 
    vplus <- 1 
    n <- 100 
    u <- runif(n, min=0, max=uplus) 
    v <- runif(n, min=0, max=vplus) 
    xi <- v/u 
    while(v < sqrt(xi)) { 
    if(v^2 <= xi) 
     return(xi) 
    } 
} 


myx <- myrnorm(1000) 
hist(myx) 

しかし、私は本当に行く方法を知りません。私はこの練習をしました。私は本当に助言に感謝します。

Answer 1

このlinkとあなたのサンプルコードの8ページに例1に続いて、私はこのソリューションを思い付いた：

ratioU <- function(nvals) 
{ 
    h_x = function(x) exp(-x) 
    # u- is b-, u+ is b+ and v+ is a in the example: 
    uminus = 0 
    uplus = 2/exp(1) 
    vplus = 1 
    X.vals <- NULL 
    i <- 0 
    repeat { 
    i <- i+1 
    u <- runif(1,0,vplus) 
    v <- runif(1,uminus,uplus) 
    X <- u/v 
    if(v^2 <= h_x(X)) { 
     tmp <- X 
    } 
    else { 
     next 
    } 
    X.vals <- c(X.vals,tmp) 
    if(length(X.vals) >= nvals) break 
    } 
    answer <- X.vals 
    answer 
} 

sol = ratioU(1000) 
par(mfrow=c(1,2)) 
hist(sol,breaks=50, main= "using ratioU",freq=F) 
hist(rexp(1000),breaks = 50, main="using rexp from R",freq=F) 
par(mfrow=c(1,1)) 

par(mfrow=c(1,2)) 
plot(density(sol)) 
plot(density(rexp(1000))) 
par(mfrow=c(1,1)) 

多くのコードを最適化することができるが、私はそれが十分このため、このように良いと思います目的。私はこれが役立つことを願っています