inverse transform cauchy dist r

Question

I'm trying to use the inverse cumulative distribution method to plot a histogram from the standard cauchy distribution and I'm getting a strange plot that doesn't look like the textbook standard cauchy. I think I have my inverse function correct (x = tan(pi*(x - 1/2))) so I would appreciate some help. Here is the r code that I have used:

n <- 10000
u <- runif(n)
c.samp <- sapply(u, function(u) tan(pi*(u - 1/2)))
hist(c.samp, breaks = 90, col = "blue",
    main = "Hist of Cauchy")

The resulting plot just doesn't look correct:

Any help is appreciated, thank you.

Answer 1

The histogram and sampling technique is correct.

Compare the results with the following (which uses the R Cauchy sampling function).

c.samp2 <- rcauchy(n)
hist(c.samp2, breaks = 90, col = "blue",
     main = "Hist of Cauchy 2")

The output here also look incorrect, but it is not.

First, you should note the x-axis is by default chosen based on the extreme values that you happen to encounter. As you probably know, the Cauchy distribution is extremely fat-tailed and very large, but rare, values are expected. When running 10000 samples from the Cauchy distribution, those relatively few single measurements squeeze the plot and do not show up on the plot because only very few observations are allocated to each bins in those extremes.

The default parameters of how hist chooses the bins are also poorly suited for distribution like the Cauchy. Try e.g.

hist(c.samp2, breaks = "FD", col = "blue",
     bins = 50,
     main = "Hist of Cauchy 2",
     xlim = c(-500, 500))

I suggest to read the help("hist") page carefully and play around with the parameters to get a good and useful histogram.

By tweaking the chosen x-axis ranges, using an y-axis probability scale, adding the theoretical distribution and a "rug", you get something more useful.

hist(c.samp, breaks = "FD", col = "blue",
     main = "Hist of Cauchy distribution",
     xlim = c(-50, 50),
     freq = FALSE)
curve(dcauchy, add = TRUE, col = "red")
rug(c.samp)

Note that using c.samp or c.samp2 now hardly changes the plot.