Inverse transform sampling : comparing bias, variance and mse for an estimator

Question

Starting from the PDF of the Pareto distribution,

\begin{equation} f_{\theta_1, \theta_2}(x) = \begin{cases} \frac{\theta_1 \theta_2^{\theta_1}}{x^{\theta_1 + 1}}, &\quad x \geq \theta_2 \\ 0, &\quad \text{otherwise} \end{cases} \end{equation}

I computed the the Gini coefficient, $$ G_{\theta_1, \theta_2} = \frac{1}{2\theta_1 - 1} $$

In particular, we can find the following estimators : MLE and MME, $$ \hat{G}_{\text{MLE}} = \frac{1}{\left( \frac{2n}{\sum_{i=1}^{n} \ln \left(\frac{X_i}{X_{(1)}} \right)} \right) - 1} $$

$$ \hat{G}_{\text{MME}} = \frac{1}{\left( \frac{2(n\bar{X} - X_{(1)})}{n(\bar{X} - X_{(1)})} \right) - 1} $$

Setting $\theta_1^0 = 3$ and $\theta_2^0 = 1$, we can generate a sample from the density $f_{\theta_1^0, \theta_2^0}$ using the inverse transform sampling. We can compute the inverse CDF,

\begin{align*} F^{-1}_{\theta_1^0, \theta_2^0}(y) &= \frac{1}{(1-y)^{1/3}} \end{align*}

Using the inverse transform sampling, I'm trying to generate $1000$ times a sample of size $20$, $40$, $60$, $80$, $100$, $150$, $200$, $300$, $400$ and $500$. For each size, we get a big sample of samples from which I would like to compute the $\textbf{bias}$, $\textbf{variance}$ and $\textbf{mean squared error}$ of the two estimators $\hat{G}_{\text{MLE}}$ and $\hat{G}_{\text{MME}}$.

Here's what I did in R,

theta_1 <- 3
theta_2 <- 1

number_of_samples <- 20

# set a seed for reproductability
set.seed(42)

# cumulative density function
cdf <- function(x) {
  (-1 / x^3)
}

# inverse of cumulative density function
inv_cdf <- function(y) {
  (1 / ((1 - y)^(1 / 3)))
}

# generate random variables vector from the inverse cdf
generate_random_variables_vector <- function(number_of_samples, inv_cdf) {
  # generate randoms numbers from the uniform distribution U(0,1)
  data_unif <- runif(number_of_samples)
  rv_vector <- inv_cdf(y = data_unif)
}

# maximum likelihood method for gini coefficient estimator
gini_mle <- function(rv_vector) {
  number_of_samples <- length(rv_vector)
  return (1 / ((2 * number_of_samples) / (sum(log(rv_vector / min(rv_vector)))) - 1))
} 

# method of moment for gini coefficient estimator
gini_mme <- function(rv_vector) {
  number_of_samples <- length(rv_vector)
  return (1 / ((2 * (number_of_samples) * mean(rv_vector) - min(rv_vector)) / (number_of_samples * (mean(rv_vector) - min(rv_vector))) - 1))
}

theoretical_gini <- function(theta_1) {
  return (1 / ((2 * theta_1) - 1))
}


n_vector <- c(20, 40, 60, 80, 100, 150, 200, 300, 400, 500)

gini_mle_biases <- numeric(10)
gini_mle_variances <- numeric(10)
gini_mle_mses <- numeric(10)

gini_mme_biases <- numeric(10)
gini_mme_variances <- numeric(10)
gini_mme_mses <- numeric(10)

i <- 0
for (n in n_vector) {
  number_of_iterations <- 1000

  gini_mle_sample <- numeric(number_of_iterations)
  gini_mme_sample <- numeric(number_of_iterations)

  for (i in 1:number_of_iterations) {
    # generate random variables
    rv_vector <- generate_random_variables_vector(number_of_samples = n, inv_cdf = inv_cdf)

    # compute gini coefficients
    gini_mle_temp <- gini_mle(rv_vector = rv_vector)
    gini_mme_temp <- gini_mme(rv_vector = rv_vector)

    gini_mle_sample[i] <- gini_mle_temp
    gini_mme_sample[i] <- gini_mme_temp
  }

  gini_mle_sample_bias <- mean(gini_mle_sample) - theoretical_gini(theta_1)
  gini_mme_sample_bias <- mean(gini_mme_sample) - theoretical_gini(theta_1)

  gini_mle_sample_variance <- sd(gini_mle_sample)^2
  gini_mme_sample_variance <- sd(gini_mme_sample)^2

  gini_mle_sample_mse <- mean((gini_mle_sample - theoretical_gini(theta_1))^2)
  gini_mme_sample_mse <- mean((gini_mme_sample - theoretical_gini(theta_1))^2)

  gini_mle_biases[i] <- mean(gini_mle_sample_bias)
  gini_mle_variances[i] <- mean(gini_mle_sample_variance)
  gini_mle_mses[i] <- mean(gini_mle_sample_mse)

  gini_mme_biases[i] <- mean(gini_mme_sample_bias)
  gini_mme_variances[i] <- mean(gini_mme_sample_variance)
  gini_mme_mses[i] <- mean(gini_mme_sample_mse)

  i <- i + 1
}

par(mfrow = c(1, 2))
plot(x = n_vector, y = gini_mle_biases, main = "", xlab = "Gini MLE biases", col = "steelblue")
plot(x = n_vector, y = gini_mme_biases, main = "", xlab = "Gini MME biases", col = "red")

par(mfrow = c(1, 2))
plot(x = n_vector, y = gini_mle_variances, main = "", xlab = "Gini MLE biases", col = "steelblue")
plot(x = n_vector, y = gini_mme_varigini_mle_variances, main = "", xlab = "Gini MME biases", col = "red")

par(mfrow = c(1, 2))
plot(x = n_vector, y = gini_mle_mses, main = "", xlab = "Gini MLE biases", col = "steelblue")
plot(x = n_vector, y = gini_mme_mle_mses, main = "", xlab = "Gini MME biases", col = "red")

I'm wondering if it's the right way to do that ?

Answer 1

Basically, what you are trying to do is evaluate the MSE : $\mathbb{E}(\hat{\theta} - \theta)^2$ ,through Monte-Carlo simulations. That is, you evaluate the expectation (an integral here) through a finite sum. Thus, you have to draw different samples to evaluate you estimator several times to form a set $\{\hat{\theta}_k\}_{k=1}^N$ from which you compute $\frac{1}{N} \sum_{k=1}^N (\hat{\theta}_k - \theta)^2 $. If that's what you do in your code (for each sample size), then you are good.