Sampling distribution of the mean of the discrete-power law distribution For a certain problem I wish to generate random integers $k$ so that their distribution follows $p_k \sim k^{-\alpha}$ for $k \geq k_{\text{min}}$, $k_{\text{min}} > 0$. I am following the procedure given in this review (page 699). Now the problem is this: I want many samples of a certain size, say size $10000$. For $\alpha = 2.2$ and $k_{\text{min}} = 2$, the theoretical value of the mean is $\langle k\rangle \approx 9.36$. Thus, when I generate my samples, and take sample averages, I expect that these averages should be close to $9.36$. However, when I plot the sampling distribution for the mean (i.e. the distribution of these sample averages), I get highly skewed distribution as shown below (total $1000$ samples were generated):

As is clear, most samples give average below the theoretical mean while some have very high values compared to the theoretical mean. One may argue that this is expected anyway because of the nature of the power-laws.

But my question is, if I want to say that my results correspond to
mean value $9.36$ would that be right, if I obtain them using each of these samples? If not, can I generate the samples so that the distribution of the sample averages would be symmetric around the theoretical mean?

I can think of the following option: For a sample of $n$ points, generate $n-1$ points from the power-law, and add $n^{\text{th}}$ point manually so that the sample average would come out right. However, I am not sure if I would be really drawing from the power-law distribution then.
Any help is highly appreciated.
 A: The distribution you are dealing with is a truncated zeta distribution, with mass function given by:
$$p_K(k) = \frac{k^{-\alpha}}{\zeta (\alpha,k_\min)}
\quad \quad \quad \text{for all integers } k \geqslant k_\min,$$
where we use the Hurwitz zeta function given (for positive integer $k_\min$) by $\zeta (\alpha,k_\min) = \sum_{k=k_\min}^\infty k^{-\alpha}$.  The mean and variance for this distribution are given respectively by:
$$\begin{align}
\mathbb{E}(K) &= \frac{\zeta (\alpha-1,k_\min)}{\zeta (\alpha,k_\min)}
\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \text{for } \alpha > 2, \\[8pt]
\mathbb{V}(K) &= \frac{\zeta (\alpha,k_\min) \zeta (\alpha-2,k_\min) - \zeta (\alpha-1,k_\min)^2}{\zeta (\alpha,k_\min)^2}
\quad \quad \quad \ \text{for } \alpha > 3. \\[6pt]
\end{align}$$
With $\alpha = 2.2$ the mean of the distribution is $\mathbb{E}(K) = \zeta(1.2,2)/\zeta(2.2,2) = 9.360199$ and its variance is infinite.  This means that the distribution is not amenable to the classical central limit theorem, but it still obeys the law of large numbers.  (It might be amenable to a generalised central limit theorem that is applicable to distributions with infinite variance.  This requires you to look at the stability of the distribution.)  Consequently, the sample mean will converge towards the true mean, but the distribution of the sample mean does not converge to a normal distribution.    One would indeed expect the distribution of the sample mean to be positively skewed, owing to the occurrence of extreme positive values under a power-law distribution.
In regard to your question, the notion that the distribution of the sample mean "corresponds" to the true expected value is not really clear, so if you say that, it does not really have a clear meaning.  What you can say is that the law of large numbers holds, so the sample mean will converge to the true mean as $n \rightarrow \infty$.

Implementation in R: For the sake of replication, I will repeat your simulation analysis to see if I get the same results you are getting.  I recommend you code your simulation so that you get a "replicable analysis" by setting the seed, etc.  The zeta distribution is contained in the VGAM package in R, which contains all the standard probability functions.  In particular, this allows us to generate values from the zeta distribution, and we can then truncate by ignoring values below the stipulated minimum.  In the code below I generate $m=1000$ samples each containing $n=10000$ data points from your distribution.
#Set parameters
kmin  <- 2;
alpha <- 2.2;
n     <- 10000;
m     <- 1000;

#Compute true mean parameter
mean.par <- VGAM::zeta(alpha-1, shift = 2)  /VGAM::zeta(alpha, shift = 2);

#Create matrix of values from truncated zeta distribution
set.seed(1);
VALUES  <- numeric(n*m);
IND     <- 0;
while (IND < n*m) {
    RAND <- VGAM::rzeta(10000, shape = alpha-1);
    RAND <- RAND[RAND >= kmin];
    RR   <- length(RAND);
    VALUES[(IND+1):(IND+RR)] <- RAND;
    IND  <- IND+RR; }
VALUES  <- VALUES[1:(n*m)];
SAMPLES <- matrix(VALUES, nrow = n, ncol = m);

#Compute sample means and plot their distribution
MEANS <- colMeans(SAMPLES);
TITLE <- paste0('Histogram of sample means \n (', m, ' samples with n = ', n, ' values)');
hist(MEANS, freq = FALSE, breaks = 150, xlim = c(0,60),
     main = TITLE, xlab = 'Sample mean');
abline(v = mean.par, col = "red", lwd = 2, lty = 2);


