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

Starting from the PDF of the Pareto distribution,

$$$$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}$$$$

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 ?

• Are you asking us to review your R code? Programming questions are off-topic on Cross Validated and are better suited to Stack Overflow but even on SO it's not appropriate to dump code and ask if the code is correct. Commented Apr 23, 2022 at 11:57

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.