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: