3
$\begingroup$

In R, suppose you generate $10^7$ random variates of the standard Cauchy distribution. Then you plot a histogram and density of this simulated data with the actual standard Cauchy density on top of it. For some reason, the charts don't match:

n = 10^7
tibble(cauchy_real = rcauchy(n)) %>%
  filter(cauchy_real >= -10, cauchy_real <= 10) %>%
  ggplot() +
  geom_density(aes(x = cauchy_real)) +
  geom_function(fun = dcauchy, colour = 'red') +
  geom_hline(yintercept = 1/pi)

R output:

enter image description here

Observations: 1. We filter the sample generated via rcauchy because it can generate some extreme outliers, and 2. We plot the horizontal bar on $y = 1/ \pi$ because, generally, in a Cauchy distribution, the peak of the distribution is at $y = 1/(\pi \gamma)$. 3. This discrepancy doesn't go away even if you increase n.

Why does this happen?

Improve this question
$\endgroup$
5
  • 7
    $\begingroup$ Hint: Please explain in more detail what you believe filter does here. Another hint: do the same thing with, say, a Normal distribution, but filter out all values more extreme than one SD from the mean. $\endgroup$
    – whuber ♦
    Commented Apr 27, 2023 at 19:10
  • $\begingroup$ Filter simply removes extreme data. In the particular case of the cauchy distribution, if you run rcauchy(10000) you get values such as -4000 which completely distort the graph. Another option is to use xlim(-5, 5) but the same thing happens. The cauchy doesn't have a closed form for SD so you cannot really compare the two. $\endgroup$
    – Sigma
    Commented Apr 27, 2023 at 19:16
  • 1
    $\begingroup$ @Sigma xlim in ggplot deletes observation outside of it's range. You could use coord_cartesian(xlim = c(-10, 10)) but geom_density doesn't cope well with extreme range you get from 10^7 Cauchy values $\endgroup$
    – Lukas Lohse
    Commented Apr 27, 2023 at 21:08
  • 3
    $\begingroup$ See what happens if you replace geom_function(fun = dcauchy, colour = 'red') with geom_function(fun = (\(x)(1/(2*atan(10) + 2*x^2*atan(10)))), colour = 'red'). The second function is the density of a truncated Cauchy distribution at $-10$ and $+10$. $\endgroup$
    – COOLSerdash
    Commented Apr 28, 2023 at 7:57
  • $\begingroup$ Your red curve should have a peak at 1/pi $\approx 0.32$ while your black curve should have a peak at 1/pi/(pcauchy(10)-pcauchy(-10)) $\approx 0.34$ and that seems to be what you see $\endgroup$
    – Henry
    Commented Apr 28, 2023 at 10:04

2 Answers 2

Reset to default
17
$\begingroup$

There is nothing surprising here since you're comparing a kernel density estimate based on a sample from a truncated Cauchy distribution with the theoretical pdf of an untruncated Cauchy.

If what you're trying to do is to plot a kernel density estimate based on a sample from an untruncated Cauchy over a range not including the whole sample and compare that to the theoretical pdf, this can be done straightforwardly in base R:

n <- 10^6
x <- rcauchy(n)
f <- density(x, from = -10, to = 10)
plot(f, col="red", main="")
curve(dcauchy, add=TRUE, lty=2)

From what I can tell (but I'm not an expert on tidyverse), there is no simple way to specify this range in the geom_density function in the tidyverse similar to how the to and from arguments in the density function is used in the above code. Adding coord_cartesian(xlim = c(-10, 10)) (as suggested by @LukasLohse) and something like n = 1e+6 in the calls to geom_density and geom_function would potentially work but would make the code painfully slow and memory intensive.

Improve this answer
$\endgroup$
5
  • 1
    $\begingroup$ This doesn't appear to address the question. To match the OP's approach you should be plotting the truncated density, as in x <- rcauchy(1e6); a <- 10; hist(x[abs(x) <= a], freq = FALSE, breaks = seq(-a, a, length.out = 100)); curve(dcauchy(x) / diff(pcauchy(c(-a, a))), add = TRUE, lwd = 2, col = "Red"); abline(h = 1/pi) Setting a to smaller values ought to make the situation abundantly clear. $\endgroup$
    – whuber ♦
    Commented Apr 27, 2023 at 22:15
  • $\begingroup$ Maybe @hadley can comment on the "proper" tidyverse way of doing this. $\endgroup$
    – Jarle Tufto
    Commented Apr 28, 2023 at 12:31
  • $\begingroup$ "Adding coord_cartesian(xlim = c(-10, 10)) (as suggested by @LukasLohse) and something like n = 1e+6 in the calls to geom_density and geom_function would work but would make the code painfully slow and memory intensive." I tried this and actually got just a weird chart that looks nothing like the cauchy or even a density. Anyway, thanks for providing a method using base R graphics. $\endgroup$
    – Sigma
    Commented Apr 28, 2023 at 12:49
  • 1
    $\begingroup$ @Sigma Yes, I should add that the (in my mind) ugly tidyverse "solution" I propose would only sometime work. You may need to increase n (the number of evaluation points of the density estimate and the theoretical pdf) to an even larger number depending on how large the difference is between the largest and smallest value in the simulated random sample. $\endgroup$
    – Jarle Tufto
    Commented Apr 28, 2023 at 12:56
  • 2
    $\begingroup$ You got the weird chart because the software chose a huge radius for the kernel and must have discretized it over huge bins. I'll bet the part of the plot you drew was composed of just a few connected line segments. But again, I think the issue raised by your question concerns what the density of a truncated distribution is and the code I offered shows how to compute the proper normalization factor. $\endgroup$
    – whuber ♦
    Commented Apr 28, 2023 at 13:13
3
$\begingroup$

If it is too difficult to scale the black line you can scale the red one and the horizontal line instead. In the following code you get overlap:

n = 10^7
tibble(cauchy_real = rcauchy(n)) %>%
  filter(cauchy_real >= -10, cauchy_real <= 10) %>%
  ggplot() +
  geom_density(aes(x = cauchy_real)) +
  geom_function(fun = function(x){dcauchy(x) / (pcauchy(10) - pcauchy(-10))}, colour = 'red') +
  geom_hline(yintercept = 1/pi / (pcauchy(10) - pcauchy(-10)))

enter image description here

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.