Best way to plot a heavy tailed distribution? Log-log seems more conventional to plot a probability distribution to look for evidence of a heavy tail. Why is this the case? For data with a heavier tail than an exponential distribution, wouldn't log-linear be best, i.e., $log(f(x))$ vs. $x$? On that plot, an exponential function would be linear.
 A: The classical concept of "heavy tails" that we are usually interested in is when the tails are sufficiently heavy to give infinite variance.  This occurs under power-law tail behaviour when one or both tails decay at a rate that is no faster than cubic decay.  (Conversely, the variance will be finite if both tails decay faster than cubic decay.)  You can find a detailed explanation of the log-log properties of a distribution with power-law tails in this related question.
If a distribution has a power-law tail of the form $f(x) \rightarrow c x^{-\omega-1}$ then it can be shown that the log-tail-probability and log-deviation from a fixed point $x_0$ are related in the tail by:
$$\ln \mathbb{P}(X>x) \rightarrow \text{const} - \omega \ln(x-x_0) \quad \quad \text{as } x \rightarrow \infty.$$
Now, if we let $x_{(1)} \geqslant \cdots \geqslant x_{(n)}$ be the ordered sample values in our data, we can estimate the log-tail-probability by $\ln \hat{\mathbb{P}}(X>x_{(i)}) = \ln(2i-1)- \ln(2n)$.  For values in the tails of the distribution we should therefore expect the values to follow the relationships:
$$\begin{matrix}
\text{Right tail } & & & & \quad \quad \quad \ln(2i-1) \approx \text{const} - \omega \ln|x_{(i)}-x_0|, \\
\text{Left tail} \quad & & & & \text{ } \ln(2(n-i)-1) \approx \text{const} - \omega \ln|x_{(i)}-x_0|. \\
\end{matrix}$$
As you can see, this gives us a log-log comparison that allows us to see whether the data exhibits power-law decay in the tails, and lets us estimate the rate of decay.  If $\omega \leqslant 2$ then there is decay no faster than cubic decay, which means the distribution has infinite variance --- i.e., it is heavy tailed.

Generating the tail-plot: You can generate the tail-plot for a set of data using the tailplot function in the utilities package in R.  This function takes in a set of data and produces the tail-plot for one or both tails.  (The function can also produce the Hill-plot and/or the DeSousa-Michailidis-plot, both of which are also useful for estimating rates of tail decay.  Below we generate some data from an exponential distribution and show the tail-plot for both tails.  We can see that the left-tail decays extremely rapidly (which is unsurprising, since it is bounded) and the right-tail also decays faster than cubic decay.  This plot provides strong empirical evidence that the distribution is not heavy tailed ---i.e., it has finite variance.
#Generate data from an exponential distribution
set.seed(1)
DATA <- rexp(1000, rate = 1)

#Generate the tail-plot for the data (both tails)
library(utilities)
tailplot(DATA)


