Using the acceptance-rejection method for simulating

768 Views Asked by At

If we have a function for which an inverse is not easily calculable, we need to use the acceptance-rejection method to simulate from it.

The problem I have is to simulate from the CDF $F_X(x)=(\frac{1}{\pi}\tan^{-1}(x))+\frac{1}{2}$.

I know that there is an algorithm, or rather set of steps to follow when using the acceptance-rejection method, but I do not at all understand it. Could someone explain it to me?

Bear in mind the question is not asking me to literally simulate from this CDF but rather to describe the algorithm which would allow one to simulate from it.

1

There are 1 best solutions below

0
On

Comment illustrated: Using three methods to simulate standard Cauchy: (a) Quantile (inverse CDF) method, (b) Ratio of two standard normals, (c) as Student's t distribution with one degree of freedom (an R function).

Below three simulations of a standard Cauchy distribution are shown, truncated to the interval $(-6, 6)$ to make useful histograms (unencumbered by very many far outliers). In each case, about 10.5% of the observations are beyond the truncation limits. The red curve is the Cauchy density, adjusted for truncation.

enter image description here

R code for simulations and histograms:

u = runif(10^6);  x1 = tan(pi*u - pi/2); X1=x1[x1 < 6 & x1 > -6]
z1 = rnorm(10^6); z2 = rnorm(10^6); x2 = z1/z2; X2=x2[x2 < 6 & x2 > -6]
x3 = rt(10^6, 1); X3=x3[x3 < 6 & x3 > -6]

K = diff(pt(c(-6,6), 1))
par(mfrow=c(3,1)) # enables three panels per plot
 hist(X1, prob=T, br=60, col="skyblue2", main="By Quantile Method")
  curve(dt(x, 1)/K, add=T, col="red")
 hist(X2, prob=T, br=60, col="skyblue2", main="As Ratio of Standard Normals")
  curve(dt(x, 1)/K, add=T, col="red")
 hist(X3, prob=T, br=60, col="skyblue2", main="As T(df=1) from R")
  curve(dt(x, 1)/K, add=T, col="red")
par(mfrow=c(1,1))

Standard Cauchy is an extremely heavy-tailed distribution. Minimums and Maximums before truncation:

range(x1); range(x2); range(x3) # in R 'range' returns min and max
[1] -1223930.7   228502.5
[1]  -79002.49 2180915.71
[1] -175184.4  330493.7

The population mean does not exist for the standard Cauchy distribution. Sample medians of un-truncated data are near 0 in all three simulations:

median(x1); median(x2); median(x3)
[1] 0.001211143
[1] -0.0009669604
[1] 0.001494456