Write a Non-Uniform RNG

Original Assignment

The distribution having density

f(x) = 1 / [π cosh(x)],     − ∞ < x < ∞
where
cosh(x) = ½ (ex + e− x)
is of no particular interest in statistics. It has the curious property that its density is proportional to its characteristic function (so the normal distribution is not the only distribution having this property).
  1. Write an R function to sample this distribution using the inverse CDF method. Since this distribution has no parameters, there will only be one argument, the sample size.
  2. Write a C function callable from R that does exactly the same thing as your R function — it produces bit for bit identical results when using the same .Random.seed (which implies that it uses the R random number generators).
  3. Provide an R batch file (a file to be run by R CMD BATCH) that calls both methods for 100 random variates and shows they produce the same results.
  4. Repeat 1, 2, and 3 where the method is rejection sampling from the double exponential distribution, using the method X = − log(U) to generate the exponentials.

The R functions 

foo <- function(x) 1 / 2 + 2 * atan(tanh(x / 2)) / pi
bar <- function(p) 2 * atanh(tan((2 * p - 1) * pi / 4))

calculate the CDF and inverse CDF, respectively, of the distribution in question. These were derived (in class) by the Mathematica code 

In[1]:= Integrate[ 1 / (Pi Cosh[x]), {x, 0, y}, Assumptions -> y > 0 ]

                      y
        2 ArcTan[Tanh[-]]
                      2
Out[1]= -----------------
               Pi

In[2]:= foo[y_] = 1 / 2 + %

                          y
            2 ArcTan[Tanh[-]]
        1                 2
Out[2]= - + -----------------
        2          Pi

In[3]:= Solve[ foo[y] == p, y ]

Solve::ifun: Inverse functions are being used by Solve, so some solutions may
     not be found; use Reduce for complete solution information.

                             (-1 + 2 p) Pi
Out[3]= {{y -> 2 ArcTanh[Tan[-------------]]}}

Comment Added Later

It occurred to me to me after some discussion in class that we did not discuss sharp tests for when random numbers are correct. A Q-Q plot is the best test. Other tests, such as chi-square and Kolmogorov-Smirnov, are not as sensitive to the tails of the distribution, although Kolmogorov-Smirnov is worth doing too.

If foo is the function calculating the CDF, then 

ks.test(x, "foo")
does the appropriate K-S test. And if bar is the function calculating the inverse CDF, then 
xx <- bar(ppoints(length(x)))
plot(xx, sort(x), xlab = "Theoretical Quantiles", ylab = "Sample Quantiles")
abline(0, 1)
does the appropriate Q-Q plot.

Comment Added Later