Estimating model parameters using other optimisation algorithms (e.g., genetic algorithm)

Question

I am trying to better understand how model parameters are estimated, so I decided to try and do this myself. Suppose we have the likelihood function of the model - I am interested in knowing if "other algorithms" (e.g. the Genetic Algorithm) can be used to estimate model parameters instead of the usual choice of the Newton-Raphson algorithm.

Using the R programming language, given some observed data, I tried to write a function that estimates the parameters (mu1, mu2, sigma1, sigma2, sigma12) of a bivariate normal distribution (i.e. fit a bivariate normal distribution to some observed data):

Suppose I have the following data (5 observations, 2 variables - to keep it simple):

set.seed(123)

x1 = rnorm(5,1,1)
x2 = rnorm(5,2,1)

data_frame = data.frame(x1, x2)

       x1      x2
1 0.4395244 3.7150650
2 0.7698225 2.4609162
3 2.5587083 0.7349388
4 1.0705084 1.3131471
5 1.1292877 1.5543380

I attempted to write the likelihood function of the normal bivariate distribution in several parts (Disclaimer: I am still learning R and this is probably a very inefficient way to write the likelihood function):

rho = sigma12/(sigma1*sigma2)
a = -1/(2 * (1 - ((sigma12 / (rho))^2))) * ((0.4395244 - mu1)/ sigma1)^2 - 2*(sigma12 / (sigma1 * sigma2)) *((0.4395244 - mu1)/sigma1) * ((3.7150650 - mu2)/sigma2)  + ((3.7150650 - mu2)/sigma2)^2
b = -1/(2 * (1 - ((sigma12 / (rho))^2))) * ((0.7698225 - mu1)/ sigma1)^2 - 2*(sigma12 / (sigma1 * sigma2)) *((0.7698225 - mu1)/sigma1) * ((2.4609162 - mu2)/sigma2)  + ((2.4609162 - mu2)/sigma2)^2
c =  -1/(2 * (1 - ((sigma12 / (rho))^2))) * ((2.5587083 - mu1)/ sigma1)^2 - 2*(sigma12 / (sigma1 * sigma2)) *((2.5587083 - mu1)/sigma1) * ((0.7349388 - mu2)/sigma2)  + ((0.7349388 - mu2)/sigma2)^2
d = -1/(2 * (1 - ((sigma12 / (rho))^2))) *  ((1.0705084 - mu1)/ sigma1)^2 - 2*(sigma12 / (sigma1 * sigma2)) *((1.0705084 - mu1)/sigma1) * ((1.3131471 - mu2)/sigma2)  + ((1.3131471 - mu2)/sigma2)^2
e = -1/(2 * (1 - ((sigma12 / (rho))^2))) * ((1.1292877 - mu1)/ sigma1)^2 - 2*(sigma12 / (sigma1 * sigma2)) *((1.1292877 - mu1)/sigma1) * ((1.5543380 - mu2)/sigma2)  + ((1.5543380 - mu2)/sigma2)^2

Next, I define the likelihood function (I call it "fitness"):

fitness <- function(mu1, mu2, sigma1, sigma2, sigma12) {

answer = answer = 1/ ((2*pi*sigma1*sigma2*sqrt((1 - rho^2)) / ((sigma1 * sigma2) * exp( a * b * c * d * e))))

}

Finally, I calculate the model parameters using the Genetic Algorithm for the optimization:

libary(GA)

GA <- ga(type = "real-valued", 
         fitness = function(x)  fitness(x[1], x[2], x[3], x[4], x[5]),
         lower = c(-5, -5, -5, -5, -5), upper = c(5, 5, 5, 5, 5), 
         popSize = 50, maxiter = 1000, run = 100)

The final model parameter estimates (and the final value of the likelihood function) can be seen below:

summary(GA)

-- Genetic Algorithm ------------------- 

GA settings: 
Type                  =  real-valued 
Population size       =  50 
Number of generations =  1000 
Elitism               =  2 
Crossover probability =  0.8 
Mutation probability  =  0.1 
Search domain = 
      x1 x2 x3 x4 x5
lower -5 -5 -5 -5 -5
upper  5  5  5  5  5

GA results: 
Iterations             = 100 
Fitness function value = 0.183913 
Solution = 
              x1       x2        x3       x4        x5
[1,] -0.06653808 1.017177 -1.410564 2.228266 -1.004855

According to the results of the optimization, the final value of the likelihood function is "0.183913". The final parameter estimates for (mu1, mu2, sigma1, sigma2, sigma12) = (-0.06653808, 1.017177, -1.410564, 2.228266, -1.004855).

Question: Mathematically speaking, can someone please tell me if what I have done makes mathematical sense? Can the parameters of a model be estimated using any optimization algorithm (e.g. the Genetic Algorithm)?

Thanks!

Notes:

• In my example, I directly optimized the likelihood function instead of the log-likelihood function.

• In my example, I generated the data from two independent normal distributions. However, I am creating a scenario where you directly observe this data and believe that it might have come from a bivariate normal distribution.

• In my example, there are only 5 observations - also, the Genetic Algorithm has a "stochastic" (random) component. Taking into consideration both of these facts, when you re-run the Genetic Algorithm on these same 5 observations - you will get different answers for your parameter estimates (e.g. sometimes the Genetic Algorithm might return several possible "sets of potential estimates") - but the corresponding likelihood value associated with these different parameter estimates are all roughly equal (0.183)

• I am not sure if confidence intervals can still be calculated on these parameter estimates when they are estimated using the Genetic Algorithm (maybe the Genetic Algorithm can be re-run many times and we can observe the "spread" of the different parameter estimates... but I am not sure if this is mathematically correct).

Reference:

• https://online.stat.psu.edu/stat505/lesson/4/4.2
• https://cran.r-project.org/web/packages/GA/vignettes/GA.html

Answer 1

In principle, you can certainly use the genetic algorithm (or any other optimisation algorithms) to try to compute the MLE. In the present case you have the benefit of having a known analytic form for the MLE:

$$\hat{\boldsymbol{\mu}} = \begin{bmatrix} \bar{x}_1 \\ \bar{x}_2 \end{bmatrix} \quad \quad \quad \quad \quad \hat{\boldsymbol{\Sigma}} = \frac{n-1}{n} \begin{bmatrix} s_1^2 & s_{1,2} \\ s_{1,2} & s_2^2 \\ \end{bmatrix},$$

where the quantities in the matrices refer to the sample mean, sample variance and sample covariance. If we apply this to your data, and use the invariance property of the MLE, we get the estimates:

#Compute the sample statistics
n       <- length(x1)
DATA    <- cbind(x1, x2)
MEAN    <- colMeans(DATA)
COV     <- cov(DATA)
COV.ADJ <- (1-1/n)*COV

#Compute the MLE
MLE <- data.frame(mu1 = MEAN[1], mu2 = MEAN[2], rho  = COV[1,2]/sqrt(COV[1,1]*COV[2,2]),
                  sigma1 = sqrt(COV.ADJ[1,1]), sigma2 = sqrt(COV.ADJ[2,2]))
row.names(MLE) <- 'MLE'

#Show the MLE
MLE

        mu1      mu2        rho sigma1   sigma2
MLE 1.19357 1.955681 -0.8219342 0.7254 1.040657

As you can see, your estimates are miles away from the true values of the MLE, computed analytically. (Also, your estimate for $\sigma_1$ is negative, which is an impossible value.) So, either you've made an error in your code, or the genetic algorithm is just shitting-the-bed in this particular case. (Note that optimisation of the likelihood instead of the log-likelihood is quite unstable.)

If you would like to explore the performance of different kinds of optimisation methods in this problem, the best way to do this would be to formulate code to apply each different optimisation method and then run these to see how accurately they approximate the MLE values (computed analytically) and how long they took to run. This will give you an idea of the accuracy and efficiency of different methods as they apply to this example, which has the benefit of having a known analytic form to compare to.

Answer 2

I am interested in knowing if "other algorithms" (e.g. the Genetic Algorithm) can be used to estimate model parameters instead of the usual choice of the Newton-Raphson algorithm.

First of all, the Newton-Raphson algorithm is not the "usual choice". There are many different optimization algorithms, we use different algorithms for different problems. For example, generalized linear models commonly are estimated using iteratively reweighted least squares algorithm.

Re-phrasing your question, one could make it more general and ask if you could use any optimization algorithm for some problem and expect it to produce exactly the same solution or equally good solution? Obviously not! We have many different optimization algorithms because they serve different purposes. For some problems, some algorithms will work better than others.

As for things like calculating confidence intervals, there is no single answer. For example, for Bernoulli distribution with mean $p$ variance is $p(1-p)$, so knowing the mean directly translates to knowing the variability. In other cases, calculating confidence itervals may be more complicated.