Alternatives for RMSE to Evaluate Goodness of Fit for Stable Distribution Parameters

Question

I am estimating the parameters (alpha, beta, gamma, delta) of a stable distribution from a list of numerical data. I used a package to generate data from one type of stable distribution, specifically the standard Levy distribution.

if (!require("stable")) {
  install.packages("stable", repos = "http://R-Forge.R-project.org")
  library(stable)
}

# Set seed for reproducibility
set.seed(123)

# Generate 100 points from a standard Lévy distribution
levy_data <- rstable(n = 100, alpha = 0.5, beta = 1, gamma = 1, 
                     delta = 0, pm = 0)

After that, I used a simple MLE method to find the four parameters.

# Perform Maximum Likelihood Estimation (MLE) to estimate parameters
levy_fit <- stable.fit(levy_data, method = 1)  
# Using method = 1 for MLE

# Print the estimated parameters
print(levy_fit)

Now since for a actual dataset, we would not have know the parameters before hand. So we need a way to measure the goodness of fit of the theoretical distribution we obtained to our data. My professor suggested using RMSE to assess the goodness of fit, similar to machine learning model evaluations. However, since there's no dependent variable and RMSE compares predicted versus actual observations, I'm questioning its suitability for this task.

I believe log-likelihood could be a better fit for assessing goodness in this context, and I am looking into other methods as well. I would appreciate any insights on:

• Why RMSE may not be suitable for probability distribution fitting.
• Alternative methods for measuring goodness of fit in stable distributions.

Answer 1

You don't get to know how good your empirical estimates are. This is why we develop estimators with desirable theoretical properties, so we have some reasonable expectation of acceptable performance when we go to the data. But we don't get to know how good our estimates are.

I think your assignment wants you to do something like I do below to calculate the RMSE of two competing estimators of the standard deviation of a Gaussian: the usual sample standard deviation vs the interquartile range.

set.seed(2024)

# Define sample size
#
N <- 100

# Define number of times to repeat draws from a distribution
#
R <- 1000

# Set standard deviation value for a Gaussian
#
set_sd <- 1

# Blank vectors to hold estimates
#
sds <- iqrs <- rep(NA, R)

# Loop R-many times
#
for (i in 1:R){
  
  # Draw from a Gaussian with a mean of 0 and the specified standard deviation
  #
  x <- rnorm(N, 0, set_sd)
  
  # Estimate the standard deviation by using the usual
  # sample standard deviation
  #
  sds[i] <- sd(x)
  
  # Estimate the standard deviation by using the interquartile range
  # This will turn out to be a bad estimator, but it is a valid statistic
  # Even constants can be estimators...weird, right? Mathematicalmonk has
  # a nice video where he examines the MSE of a constant estimator:
  # https://www.youtube.com/watch?v=Ow7xMTF2BfQ
  #
  iqrs[i] <- IQR(x)
}

# Since we know the standard deviation in this simulation, we can estimate
# the RMSE of each estimator
#
rmse_sd <- sqrt(mean((sds - set_sd)^2))
rmse_iqr <- sqrt(mean((iqrs - set_sd)^2))

# The RMSE of the usual sample standard deviation is much lower than that of
# the interquartile range.
#
rmse_sd  # I get 0.07069313
rmse_iqr # I get 0.3630151

In your case, you would estimate the parameters of your stable distribution instead of the standard deviation of the Gaussian that I did, but the idea is the similar. If you want to consider the RMSE of your entire vector estimator, perhaps you can use the calculation described here.