Se trabajará con la función de la distribución normal:

\(f(x_{j})=\frac{1}{\sqrt{2\pi\sigma_{0}^{2}}}exp(-\frac{1}{2}\frac{(x_{j}-\mu_{0})^{2}}{\sigma_{0}^{2}} )\)

Se desarrolla la productoria:

\(L(\mu,\sigma^{2};x_{1},...,x_{n})=\prod_{j=1}^{n}f(x_{j},\mu,\sigma^{2}) =\prod_{j=1}^{n}\frac{1}{\sqrt{2\pi\sigma^{2}}}exp(-\frac{1}{2}\frac{(x_{j}-\mu)^{2}}{\sigma^{2}})=(2\pi\sigma^{2})^{-n/2}exp(-\frac{1}{2\sigma^{2}}\sum_{j=1}^{n}(x_{j}-\mu)^{2})\)

Así se obtiene la siguiente función L(likelihood function):

\(L(\mu,\sigma^{2};x_{1},...,x_{n})=(2\pi\sigma^{2})^{-n/2}exp(-\frac{1}{2\sigma^{2}}\sum_{j=1}^{n}(x_{j}-\mu)^{2})\)

Luego aplicando logaritmo natural a ambos lados:

\(l(\mu,\sigma^{2};x_{1},...,x_{n})= ln(L(\mu,\sigma^{2};x_{1},...,x_{n})) =-\frac{n}{2}ln(2\pi)-\frac{n}{2}ln(\sigma^{2})-\frac{1}{2\sigma ^{2}}\sum_{j=1}^{n}(x_{j}-\mu)^{2}\)

De esta manera se obtiene la función l(log.likelihood function) [1]:

\(l(\mu,\sigma^{2};x_{1},...,x_{n})=-\frac{n}{2}ln(2\pi)-\frac{n}{2}ln(\sigma^{2})-\frac{1}{2\sigma ^{2}}\sum_{j=1}^{n}(x_{j}-\mu)^{2}\)

Se utilizará esta función para estimar los valores de media y desviacion estandar asociados a los datos, los cuales tendran una distribución normal.

PASO 1: Generar los datos

x=rnorm(n=1000,mean=10, sd=2)  #Se genera una lista de valores que siguen una distribucion normal utilizando la funcion rnorm. 
#Se generaran 1000 valores aleatorios que seguiran una distribucion normal de media 10 y desviacion estandar 2

Los valores de media y desviación estandar con los cuales fueron generados los datos son 10 y 2 respectivamente, al optimizar la funcion l obtenida anteriormente se deberán obtener valores cercanos.

PASO 2: Derivar la función de verosimilitud

A continuacion se define una funcion en R la cual reciba un vector de valores y otro vector con la media y la desviacion estandar, en ese orden.

mle=function(ValoresX,ParametrosDist) {   #Se define la funcion de maxima verosimilitud, recibe los valores que siguen la distribucion normal y un vector de la forma c(media,distribucion), esta funcion sera maximizada en el paso 3 para estimar valores de media y desviacion estandar.
  n=length(ValoresX)
  media= ParametrosDist[1]
  desviacion=ParametrosDist[2]
 resultado = -((n/2)*log(2*pi))-((n/2)*log(desviacion^2))-((1/(2*desviacion^2)*( sum((ValoresX-media)^2) ))) #funcion l(log.likelihood function) obtenida en los antecedentes.
  return(-resultado)  #el resultado se entrega en negativo, ya que al optimizar con optim se necesita que maximice (por defecto optim niminiza)
}

PASO 3: Estimación por máxima verosimilitud

Se utiliza optim para maximizar a la funcion de verosimilitud (mle) definida en el paso 2, es importante mencionar que dicha funcion retorna un resultado negativo, ya que por defecto optim busca minimizar y deseamos que maximice.

Optimizacion = optim(fn=mle,                   # Función de verosimilitud
                       par=c(1,1),                    # Estimación inicial
                       lower = c(-Inf, -Inf),         # Límite inferior de los parámetros
                       upper = c(Inf, Inf),           # Límite superior de los parámetros
                       hessian=TRUE,                  # Devuelve el Hessiano 
                       method = "L-BFGS-B",
                       # Entradas personalizadas                     
                       ValoresX = x)

# Examinar estimaciones
MLE_par <- Optimizacion$par
EstimacionMedia=MLE_par[1]
EstimacionSd=MLE_par[2]

print(paste("La estimacion de la media utilizando Maxima Verosimilitud es: ", EstimacionMedia))
## [1] "La estimacion de la media utilizando Maxima Verosimilitud es:  10.0030648656853"
print(paste("La estimacion de la desviacion estandar utilizando Maxima Verosimilitud es: ", EstimacionSd))
## [1] "La estimacion de la desviacion estandar utilizando Maxima Verosimilitud es:  2.05171257901476"

De esta manera, se puede evidenciar que el valor estimado usando maxima verosimilitud es muy cercano a los valores de media y desviacion estandar reales, en conclusion, utilizando maxima verosimilitud se pudo estimar correctamente los datos de media y desviacion estandar de los datos. En las pruebas realizadas en este momento la estimacion de la media fue 10.0376811426738 y la desviacion estandar 2.059348981202, lo cual es muy cercano a los valores con los cuales fueron generados los datos usando rnormal (media 10 y desviacion estandar 2), es importante mencionar que rnorm genera valores aleatorios (con la media y desviacion entregadas), por lo tanto estos resultados pueden variar al ejecutarlos, sin embargo deberan ser resultados cercanos a media 10 y desviacion2

PASO 4: Graficar estimaciones

Se puede evidenciar que el maximo esta cercano del valor real (10).

mu_vals <- seq(0,20, by =1)
#Se realizan varias evaluaciones de la funcion mle dejando desviacion estandar constante 2.
Lista= list()
for (i in mu_vals) {
  Resultado = - mle(x,c(mu_vals[i],2))
  Lista = c(Lista, Resultado)
}

plot(mu_vals, Lista, type='o', main = 'Función de Verosimilitud considerando desviacion estandar 2', ylab = 'log.likelihood function', col = "blue")

Se puede evidenciar que el maximo esta cercano del valor real (2).

sd_vals <- seq(0,20, by =1)
#Se realizan varias evaluaciones de la funcion mle dejando media constante 10.
Lista2= list()
for (i in sd_vals) {
  Resultado2 = - mle(x,c(10,sd_vals[i])) #
  Lista2 = c(Lista2, Resultado2)
}

plot(sd_vals, Lista2, type='o', main = 'Función de Verosimilitud considerando media 10', ylab = 'log.likelihood function', col = "blue")

Referencias utilizadas:

1. Normal distribution - Maximum likelihood estimation. (s. f.). https://www.statlect.com/fundamentals-of-statistics/normal-distribution-maximum-likelihood https://www.statlect.com/fundamentals-of-statistics/normal-distribution-maximum-likelihood

2. Villalobos, M. (2022). Estimación puntual de parámetros y distribuciones muestrales (Parte 2). Uvirtual usach.cl https://uvirtual.usach.cl/moodle/course/view.php?id=24173§ion=24