The literature on the mathematical theory underlying Zipf's law is quite vast, and includes a large number of underlying theoretical models in which the law emerges. Zipf's law is related to power laws through the fact that it asserts a power-law relationship for the rank versus frequency of the objects under analysis, so there is also a substantial literature examining the connections between the Zipf distribution and power-law behaviours in the Pareto distribution. The statistical literature on this topic is quite vast, but you can find a good introductory exposition on this field in Mitzenmacher (2003). As you will see from that reference, there are a number of modelling approaches that lead to the behaviour set out in Zipf's law.
For natural language and vocabulary analysis, the most prominent modelling approach is an information-theoretic derivation akin to the work of Mandelbrot (1953). This paper uses information-optimisation to derive a slightly generalised form for Zipf's law; this model has had a large impact in information theory and has led to a range of later models. The approach used by Mandelbrot leads to a slighly generalised form of the Zipf distribution over the support $1,...,N$, defined by the proportionality relationship:
$$f(k|s,c,N) \propto \frac{1}{(k+c)^s} \cdot \mathbb{I}(k \in \{ 1,...,N \}),$$
with parameters $c \geqslant 0$ and $s > 0$. This distributional relationship is often exhibited on a log-log plot via the fact that the distribution satisfies:
$$\begin{align}
\log f(k|s,c,N)
&= \text{const} - s \log(k+c) \\[6pt]
&= \text{const} - s \log(k) - s \log \Big( 1+\frac{c}{k} \Big). \\[6pt]
\end{align}$$
In the special case where $c=0$ we see that the rank-frequency relationship will appear as a negative linear relationship on a log-log plot. For $c>0$ the relationship will appear nonlinear, but will become linear for $k \gg c$ (i.e., it will be close to linear except when $k$ is relatively low).
A useful starting point to investigate Zipf's law in empirical data is to plot the rank versus frequency on a log-log plot to see if it appears to roughly follow the above form. You can do this using the zipfplot function in the utilities package in R (example shown below). Inspection of the Zipf plot should give you a reasonable idea of whether or not your data follow the expected form under that distribution.
Once you have inspected your data on a Zipf plot, you can obtain the MLEs for the parameters $c$ and $s$ and use this to superimpose the estimated Zipf distribution onto the log-log plot, to see how closely the data follows the closest version of this distribution. You can also use goodness of fit tests to see if the variation of the data from the theoretical distribution is sufficient to falsify the assumed distributional form.
Now, if your data does depart from Zipf's law (in the generalised sense shown here) then that means you will need to investigate broader distributional forms that accomodate your data. It is a bad idea to try to "justify" Zipf's law against the evidence provided by your data --- you should allow your data to lead the analysis and seek models and distributional forms that are plausible when compared with your data. If your data does not fit the family of Zipf distributions, ideally you would broaden your analysis by examining distributional families that arise from some plausible simple change to the underlying information-theoretic models. Ideally you will end up with a distributional form that has a solid information-theoretic foundation and also fits your data well.
An example of the Zipf Plot: If you have an underlying vector of the data values from a discrete distribution, you can generate a Zipf plot of the data using the zipfplot function in the utilties package. The function automatically computes the ranks and frequencies of the outcomes in the data, so you enter the data into the function in its raw form. Here we show a simple example using data from a binomial distribution. As you can see from the plot, the distribution does not follow Zipf's law.
#Generate some mock data from a discrete distribution
set.seed(1)
XX <- rbinom(10000, size = 120, prob = 0.3)
#Plot the Zipf plot for the data
library(utilities)
zipfplot(XX)
