MÁXIMA VEROSIMILITUD

Función de Verosimilitud Recordemos que si \(x_1,x_2,…,x_n\) es una m.a.s. extraída de una población \(X\) con distribución de probabilidad dada por \(f(x|\theta)\), siendo \(\theta\) un (o más) parárametro poblacional deconocido, entonces la función de verosimilitud de la muestra se define mediante la expresión:

Ahora bien, como cada \(x_i\) es una realización de la v.a. \(X_i\) y estas son independientes e idénticamente distribuidas (iid), la función de verosimilitud puede escribirse como sigue:

\[\begin{equation} L(\theta)=\prod\limits_{i=1}^n f(x_i | \theta) \end{equation}\]

Estimación por el Método de Máxima Verosimilitud

La teoría nos indica que para encontrar el estimador máximo verosímil (EMV) para el parámetro poblacional desconocido \(\theta\), debemos resolver el siguiente problema de maximización:

\[\begin{equation} \max_{\widehat{\theta}}\; L(\theta)\quad\text{o bien} \quad \max_{\widehat{\theta}}\; \ln L(\theta) \end{equation}\]

Sin embargo, computacionalmente hablando, el problema que realmente se resuelve es minimizar el negativo de la función de log-verosimilitud:

\[\begin{align} -\ln L(\theta)&= -\ln \prod\limits_{i=1}^n f(x_i | \theta)\\ &=-\sum\limits_{i=1}^n \ln f(x_i|\theta) \end{align}\]

Lo anterior es posible gracias a que \(maxZ=min(−Z)\), siendo Z una función objetivo cualquiera.

install.packages("bbmle")

## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.1'
## (as 'lib' is unspecified)

install.packages("stats4")

## Installing package into '/cloud/lib/x86_64-pc-linux-gnu-library/4.1'
## (as 'lib' is unspecified)

## Warning: package 'stats4' is a base package, and should not be updated

library(stats4) # para la función mle
library(bbmle) # para la función mle2

Función mle del paquete stats4

Comenzaremos generando 50 datos aleatorios según una distribución normal con media 10 y desviación estándar 2.

set.seed(123)
x = rnorm(n = 100,mean = 10,sd = 2)
x

##   [1]  8.879049  9.539645 13.117417 10.141017 10.258575 13.430130 10.921832
##   [8]  7.469878  8.626294  9.108676 12.448164 10.719628 10.801543 10.221365
##  [15]  8.888318 13.573826 10.995701  6.066766 11.402712  9.054417  7.864353
##  [22]  9.564050  7.947991  8.542218  8.749921  6.626613 11.675574 10.306746
##  [29]  7.723726 12.507630 10.852928  9.409857 11.790251 11.756267 11.643162
##  [36] 11.377281 11.107835  9.876177  9.388075  9.239058  8.610586  9.584165
##  [43]  7.469207 14.337912 12.415924  7.753783  9.194230  9.066689 11.559930
##  [50]  9.833262 10.506637  9.942906  9.914259 12.737205  9.548458 13.032941
##  [57]  6.902494 11.169227 10.247708 10.431883 10.759279  8.995353  9.333585
##  [64]  7.962849  7.856418 10.607057 10.896420 10.106008 11.844535 14.100169
##  [71]  9.017938  5.381662 12.011477  8.581598  8.623983 12.051143  9.430454
##  [78]  7.558565 10.362607  9.722217 10.011528 10.770561  9.258680 11.288753
##  [85]  9.559027 10.663564 12.193678 10.870363  9.348137 12.297615 11.987008
##  [92] 11.096794 10.477463  8.744188 12.721305  8.799481 14.374666 13.065221
##  [99]  9.528599  7.947158

hist(x, main ="HISTOGRAMA MU=10",xlab = "VALORES DE X", ylab = "FRECUENCIAS", col=5)

library(ggplot2)
qplot(x,geom="histogram",bins=8, color = "red")

qplot(x,geom="density",fill = "red")+
 ggtitle("FUNCION DE DENSIDAD")+
  labs(x="valores de x", y="frecuencia")+
  theme_minimal()

Ahora construimos el negativo de la función de log-verosimilitud. Esta función, llamada NegLogLik en este caso, depende de los parámetros poblacionales desconocidos media y desviación estándar, mu y sigma respectivamente:

NegLogLik = function(mu,sigma){-sum(dnorm(x,mu,sigma,log = TRUE))}

Luego usamos la función mle del paquete stats4 para realizar la estimación por máxima verosimilitud:

EMV1 = mle(NegLogLik, start = list(mu=10, sigma=5))

## Warning in dnorm(x, mu, sigma, log = TRUE): NaNs produced

summary(EMV1)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = NegLogLik, start = list(mu = 10, sigma = 5))
## 
## Coefficients:
##        Estimate Std. Error
## mu    10.180812  0.1816481
## sigma  1.816481  0.1284445
## 
## -2 log L: 403.1679

Por otro lado, en el ejemplo anterior se utilizaron datos generados aleatoriamente con ayuda de las distribuciones de probabilidad integradas en R.

Si solo contamos con datos sueltos como por ejemplo 9, 10, 5, 3, 6, 8 y 10, los cuales suponemos siguen una distribución normal, el proceso para realizar la estimación por máxima verosimilitud es similar al anterior.

Lo primero que hacemos es construir un vector y en el cual guardamos los datos muestrales:

y = c(9,10,5,3,6,8,10)
y

## [1]  9 10  5  3  6  8 10

length(y)

## [1] 7

Construimos el negativo de la función de log-verosimilitud respectiva:

NegLogLik2 = function(mu,sigma){-sum(dnorm(y,mu,sigma,log = TRUE))}

Realizamos la estimación por máxima verosimilitud:

EMV2 = mle(NegLogLik2, start = list(mu=10, sigma=5))

## Warning in dnorm(y, mu, sigma, log = TRUE): NaNs produced

## Warning in dnorm(y, mu, sigma, log = TRUE): NaNs produced

Y finalmente visualizamos los resultados obtenidos en la estimación:

summary(EMV2)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = NegLogLik2, start = list(mu = 10, sigma = 5))
## 
## Coefficients:
##       Estimate Std. Error
## mu    7.285717  0.9414340
## sigma 2.490800  0.6656942
## 
## -2 log L: 32.64159

## Función mle2 del paquete bbmle

La función mle2 del paquete bbmle en su forma más básica es prácticamente igual a la función mle del paquete stats4. La principal diferencia radica en que debemos indicar en el argumento data de donde se tomarán los datos para realizar los cálculos:

EMV3 = mle2(NegLogLik,start = list(mu=10,sigma=5), data = list(x))

## Warning in dnorm(x, mu, sigma, log = TRUE): NaNs produced

Visualizamos los datos y vemos que los resultados son prácticamente iguales a los obtenidos anteriormente. No obstante, la función mle2 entrega más información que la función mle, particularmente, el término Pr(z), el cual será discutido más adelante.

summary(EMV3)

## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = NegLogLik, start = list(mu = 10, sigma = 5), 
##     data = list(x))
## 
## Coefficients:
##       Estimate Std. Error z value     Pr(z)    
## mu    10.18081    0.18165  56.047 < 2.2e-16 ***
## sigma  1.81648    0.12844  14.142 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 403.1679

La principal ventaja de la función mle2 es que con ella se puede prescindir de construir la función de verosimilitud, ya que es posible pedirle a la función que la calcule internamente usando el símbolo ~.

A modo de ejemplo, generemos una secuencia de 30 números según una distribución normal con media igual a 20 y desviación típica igual a 5:

set.seed(12345)
x1=rnorm(30, mean=20, sd=5)
EMV4 = mle2(x1~dnorm(mu,sigma),start = list(mu=15,sigma=4), data=data.frame(x1))

## Warning in calc_mle2_function(minuslogl, parameters, start = start, parnames =
## parnames, : using dnorm() with sd implicitly set to 1 is rarely sensible

summary(EMV4)

## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = x1 ~ dnorm(mu, sigma), start = list(mu = 15, 
##     sigma = 4), data = data.frame(x1))
## 
## Coefficients:
##       Estimate Std. Error z value     Pr(z)    
## mu    20.39395    0.84206 24.2190 < 2.2e-16 ***
## sigma  4.61218    0.59541  7.7462 9.467e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 176.8597

data()

Notemos que la notación x~dnorm(mu,sigma) significa que \(x\) es en realidad el primer argumento de la función dnorm(). De igual forma, es importante destacar que al usar esta metodología, los datos que se utilizarán para realizar los cómputos se deben ingresar como un data.frame().

Ejercicios

Use R para implementar y resolver los siguientes problemas:

Genere 50 datos aleatorios que sigan una distribución de Poisson de parámetro \(\lambda=5\) y utilice estos datos para estimar el verdadero valor de \(\lambda\). Para hacer esto, utilice por separado los siguientes métodos: mle, mle2 y mle2 sin escribir explícitamente la función de log-verosimilitud.

set.seed(123)
z = rpois(50,5)
z

##  [1]  4  7  4  8  9  2  5  8  5  5  9  5  6  5  2  8  3  2  4  9  8  6  6 11  6
## [26]  6  5  5  4  3  9  8  6  7  1  5  6  3  4  3  3  4  4  4  3  3  3  5  4  7

barplot(z, main ="BARRA",xlab = "VALORES DE X", ylab = "FRECUENCIAS", col=10)

library(ggplot2)
qplot(z,geom="boxplot",bins=8, color = "blue")

## Warning: Ignoring unknown parameters: bins

qplot(z,geom="density",fill = "red")+
 ggtitle("FUNCION DE DENSIDAD")+
  labs(x="valores de x", y="frecuencia")+
  theme_minimal()

NegLogLik = function(lambda){-sum(dpois(z,lambda,log = TRUE))}

EMV1 = mle(NegLogLik, start = list(lambda=5))

summary(EMV1)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = NegLogLik, start = list(lambda = 5))
## 
## Coefficients:
##        Estimate Std. Error
## lambda 5.240458  0.3237565
## 
## -2 log L: 219.0201

NegLogLik2 = function(lambda){-sum(dpois(z,lambda,log = TRUE))}

EMV2 = mle(NegLogLik2, start = list(lambda=5))

summary(EMV2)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = NegLogLik2, start = list(lambda = 5))
## 
## Coefficients:
##        Estimate Std. Error
## lambda 5.240458  0.3237565
## 
## -2 log L: 219.0201

EMV3 = mle2(NegLogLik,start = list(lambda=5), data = list(z))

summary(EMV3)

## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = NegLogLik, start = list(lambda = 5), data = list(z))
## 
## Coefficients:
##        Estimate Std. Error z value     Pr(z)    
## lambda  5.24046    0.32376  16.186 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 219.0201

Haga lo mismo que el ejercicio anterior pero para los parámetros de las siguientes distribuciones de probabilidad: Weibull, gamma, chi-cuadrado y exponencial. Puede elegir los parámetros poblacionales a su antojo.
Weibull

set.seed(123)
a = rweibull(50,2,1)
a

##  [1] 1.11636142 0.48771922 0.94556679 0.35271854 0.24774668 1.75749878
##  [7] 0.79903644 0.33737142 0.77151233 0.88538990 0.21006200 0.88945251
## [13] 0.62389220 0.74666562 1.50789847 0.32489231 1.18408917 1.78007562
## [19] 1.05592775 0.21578650 0.34212801 0.60581268 0.66745451 0.07580702
## [25] 0.64964842 0.58699422 0.78018246 0.72155173 1.11390126 1.38439513
## [31] 0.19410488 0.32063886 0.60831083 0.47835697 1.92469548 0.85940181
## [37] 0.52580018 1.23717018 1.07010974 1.20939332 1.39510217 0.93838719
## [43] 0.93944420 0.99868792 1.37147842 1.40523220 1.20688462 0.87387083
## [49] 1.15080921 0.39160184

hist(a, main ="histograma",xlab = "VALORES DE X", ylab = "FRECUENCIAS", col=12)

library(ggplot2)
qplot(a,geom="histogram",bins=8, color = "red")

qplot(a,geom="density",fill = "red")+
 ggtitle("FUNCION DE DENSIDAD")+
  labs(x="valores de x", y="frecuencia")+
  theme_minimal()

NegLogLik = function(shape,scale){-sum(dweibull(a,shape,scale,log = TRUE))}

EMV1 = mle(NegLogLik, start = list(shape=2, scale=1))

## Warning in dweibull(a, shape, scale, log = TRUE): NaNs produced

## Warning in dweibull(a, shape, scale, log = TRUE): NaNs produced

summary(EMV1)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = NegLogLik, start = list(shape = 2, scale = 1))
## 
## Coefficients:
##       Estimate Std. Error
## shape 1.989822 0.22515816
## scale 0.953204 0.07121771
## 
## -2 log L: 55.99141

NegLogLik2 = function(shape,scale){-sum(dweibull(a,shape,scale,log = TRUE))}

EMV2 = mle(NegLogLik2, start = list(shape=2, scale=1))

## Warning in dweibull(a, shape, scale, log = TRUE): NaNs produced

## Warning in dweibull(a, shape, scale, log = TRUE): NaNs produced

summary(EMV2)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = NegLogLik2, start = list(shape = 2, scale = 1))
## 
## Coefficients:
##       Estimate Std. Error
## shape 1.989822 0.22515816
## scale 0.953204 0.07121771
## 
## -2 log L: 55.99141

EMV3 = mle2(NegLogLik,start = list(shape=2,scale=1), data = list(a))

## Warning in dweibull(a, shape, scale, log = TRUE): NaNs produced

## Warning in dweibull(a, shape, scale, log = TRUE): NaNs produced

summary(EMV3)

## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = NegLogLik, start = list(shape = 2, scale = 1), 
##     data = list(a))
## 
## Coefficients:
##       Estimate Std. Error z value     Pr(z)    
## shape 1.989822   0.225158  8.8374 < 2.2e-16 ***
## scale 0.953204   0.071218 13.3843 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 55.99141

-GAMMA

set.seed(123)
w = rgamma(n = 100,shape = 2,rate = 1)
w

##   [1] 0.8920936 3.3118474 0.1443750 1.6625233 4.3358790 2.1176157 0.3507177
##   [8] 0.1304001 3.3737820 1.9730466 2.0309972 1.6386208 0.8964763 3.5048576
##  [15] 2.6814547 1.4325741 0.9768341 1.2472694 0.5065776 1.3972026 0.1454646
##  [22] 0.4121286 0.4299115 1.1336216 1.1603795 1.3456917 2.6749761 2.4619650
##  [29] 2.2551140 1.4251322 1.8714764 0.3775966 1.2561618 1.3059447 3.3442383
##  [36] 0.4398217 2.0350751 0.6107837 1.3996318 1.0174353 1.4479541 1.8571056
##  [43] 0.4903600 0.2028318 0.4757009 1.7761310 2.0009930 0.9478641 5.2886410
##  [50] 0.3108821 1.8947776 2.0991656 1.5656190 2.8421867 5.0615425 0.9588900
##  [57] 0.5996164 0.7571514 2.3543424 1.4924183 0.3774702 0.6945218 1.5070680
##  [64] 2.0089804 1.0803832 0.7064538 1.9338681 3.1441119 2.0803320 1.1273748
##  [71] 1.7634013 0.9609556 1.8066337 0.8295418 3.5002587 0.9859827 3.9642808
##  [78] 1.2252159 1.7333934 4.4694851 1.5465069 3.2864430 0.8939267 1.4453594
##  [85] 2.1047045 2.2146580 1.2312427 0.8781028 2.1536899 1.1921919 2.2035873
##  [92] 1.8915094 1.6322182 0.8179246 2.9092421 2.7280885 2.0567668 1.9452526
##  [99] 2.3193992 4.6675609

hist(w, main ="histograma",xlab = "VALORES DE X", ylab = "FRECUENCIAS", col=12)

library(ggplot2)
qplot(w,geom="histogram",bins=8, color = "red")

qplot(w,geom="density",color =4)+
 ggtitle("FUNCION DE DENSIDAD")+
  labs(x="valores de x", y="frecuencia")+
  theme_minimal()

NegLogLik = function(shape,rate){-sum(dgamma(w,shape,rate,log = TRUE))}

EMV1 = mle(NegLogLik, start = list(shape=2, rate=1))

## Warning in dgamma(w, shape, rate, log = TRUE): NaNs produced

## Warning in dgamma(w, shape, rate, log = TRUE): NaNs produced

summary(EMV1)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = NegLogLik, start = list(shape = 2, rate = 1))
## 
## Coefficients:
##       Estimate Std. Error
## shape 2.193234  0.2896872
## rate  1.274026  0.1889907
## 
## -2 log L: 279.8997

NegLogLik2 = function(shape,rate){-sum(dgamma(w,shape,rate,log = TRUE))}

EMV2 = mle(NegLogLik2, start = list(shape=2, rate=1))

## Warning in dgamma(w, shape, rate, log = TRUE): NaNs produced

## Warning in dgamma(w, shape, rate, log = TRUE): NaNs produced

summary(EMV2)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = NegLogLik2, start = list(shape = 2, rate = 1))
## 
## Coefficients:
##       Estimate Std. Error
## shape 2.193234  0.2896872
## rate  1.274026  0.1889907
## 
## -2 log L: 279.8997

EMV3 = mle2(NegLogLik,start = list(shape=2,rate=1), data = list(w))

## Warning in dgamma(w, shape, rate, log = TRUE): NaNs produced

## Warning in dgamma(w, shape, rate, log = TRUE): NaNs produced

summary(EMV3)

## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = NegLogLik, start = list(shape = 2, rate = 1), 
##     data = list(w))
## 
## Coefficients:
##       Estimate Std. Error z value     Pr(z)    
## shape  2.19323    0.28969  7.5710 3.703e-14 ***
## rate   1.27403    0.18899  6.7412 1.571e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 279.8997

chi-cuadrado

set.seed(123)
v = rchisq(n=100,df=4)
v

##   [1]  1.7841871  6.6236948  0.2887499  3.3250466  8.6717580  4.2352314
##   [7]  0.7014354  0.2608002  6.7475639  3.9460933  4.0619944  3.2772415
##  [13]  1.7929525  7.0097153  5.3629094  2.8651482  1.9536682  2.4945389
##  [19]  1.0131552  2.7944052  0.2909292  0.8242572  0.8598231  2.2672431
##  [25]  2.3207590  2.6913834  5.3499522  4.9239299  4.5102280  2.8502644
##  [31]  3.7429528  0.7551932  2.5123236  2.6118895  6.6884766  0.8796435
##  [37]  4.0701503  1.2215674  2.7992635  2.0348706  2.8959082  3.7142112
##  [43]  0.9807201  0.4056635  0.9514019  3.5522620  4.0019861  1.8957283
##  [49] 10.5772819  0.6217643  3.7895551  4.1983313  3.1312380  5.6843733
##  [55] 10.1230850  1.9177800  1.1992327  1.5143029  4.7086847  2.9848366
##  [61]  0.7549404  1.3890436  3.0141359  4.0179609  2.1607665  1.4129076
##  [67]  3.8677362  6.2882238  4.1606641  2.2547496  3.5268026  1.9219113
##  [73]  3.6132674  1.6590835  7.0005174  1.9719654  7.9285617  2.4504317
##  [79]  3.4667868  8.9389702  3.0930137  6.5728860  1.7878534  2.8907188
##  [85]  4.2094090  4.4293160  2.4624854  1.7562056  4.3073799  2.3843838
##  [91]  4.4071745  3.7830187  3.2644365  1.6358493  5.8184842  5.4561769
##  [97]  4.1135336  3.8905052  4.6387985  9.3351219

hist(v, main ="histograma",xlab = "VALORES DE X", ylab = "FRECUENCIAS", col=19)

library(ggplot2)
qplot(v,geom="histogram",bins=8, color =2)

qplot(v,geom="density",color ="red")+
 ggtitle("FUNCION DE DENSIDAD")+
  labs(x="valores de x", y="frecuencia")+
  theme_dark()

NegLogLik = function(df){-sum(dchisq(v,df,log = TRUE))}

EMV1 = mle(NegLogLik, start = list(df=4))

## Warning in dchisq(v, df, log = TRUE): NaNs produced

summary(EMV1)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = NegLogLik, start = list(df = 4))
## 
## Coefficients:
##    Estimate Std. Error
## df 3.636184  0.2344659
## 
## -2 log L: 420.9726

NegLogLik2 = function(df){-sum(dchisq(v,df,log = TRUE))}

EMV2 = mle(NegLogLik2, start = list(df=4))

## Warning in dchisq(v, df, log = TRUE): NaNs produced

summary(EMV2)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = NegLogLik2, start = list(df = 4))
## 
## Coefficients:
##    Estimate Std. Error
## df 3.636184  0.2344659
## 
## -2 log L: 420.9726

EMV3 = mle2(NegLogLik,start = list(df=4), data = list(v))

## Warning in dchisq(v, df, log = TRUE): NaNs produced

summary(EMV3)

## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = NegLogLik, start = list(df = 4), data = list(v))
## 
## Coefficients:
##    Estimate Std. Error z value     Pr(z)    
## df  3.63618    0.23447  15.508 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 420.9726

Exponencial

set.seed(123)
d = rexp(n=100, rate=1)
d

##   [1] 0.843457261 0.576610271 1.329054868 0.031577359 0.056210976 0.316501216
##   [7] 0.314227292 0.145266804 2.726236464 0.029153447 1.004830058 0.480214728
##  [13] 0.281013628 0.377117831 0.188284041 0.849786130 1.563203540 0.478760416
##  [19] 0.590934835 4.041011711 0.843149731 0.965871211 1.485275794 1.348044486
##  [25] 1.168528984 1.605852343 1.496742869 1.570652547 0.031767744 0.597849691
##  [31] 2.167839745 0.506615729 0.259557817 2.596892116 1.229025732 0.790681759
##  [37] 0.629280078 1.254641003 0.588684642 1.129290034 0.420364803 7.211007576
##  [43] 0.845721965 0.225542007 1.100338818 2.248305692 1.363734299 0.576391668
##  [49] 2.725275850 1.312163043 0.090591350 0.306203850 1.067213070 0.313516256
##  [55] 0.974640151 1.887823315 0.564588603 2.576961329 1.047695748 1.024441342
##  [61] 1.027869616 0.284668461 1.563051889 0.042088293 0.098630890 0.098569312
##  [67] 0.280385161 0.295786959 0.972429649 0.924027224 1.642404527 1.619912768
##  [73] 2.536145540 1.521549623 0.380014203 0.238512967 0.466487424 0.042271452
##  [79] 0.319676899 0.643610920 0.572563103 0.216055503 4.498673398 1.857090553
##  [85] 0.685458638 1.439452626 1.731153984 1.244783300 1.463300560 1.537393979
##  [91] 0.004599127 1.108765468 0.299970318 1.192003007 1.114928704 0.067375891
##  [97] 0.480668721 1.570454339 0.259946107 1.856922287

hist(d, main ="histograma",xlab = "VALORES DE X", ylab = "FRECUENCIAS", col=21)

library(ggplot2)
qplot(d,geom="histogram",bins=8, color =2)

qplot(d,geom="density",color ="red")+
 ggtitle("FUNCION DE DENSIDAD")+
  labs(x="valores de x", y="frecuencia")+
  theme_dark()

NegLogLik = function(rate){-sum(dexp(d,rate,log = TRUE))}

EMV1 = mle(NegLogLik, start = list(rate=1))

## Warning in dexp(d, rate, log = TRUE): NaNs produced

summary(EMV1)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = NegLogLik, start = list(rate = 1))
## 
## Coefficients:
##       Estimate Std. Error
## rate 0.9562808 0.09562797
## 
## -2 log L: 208.9409

NegLogLik2 = function(rate){-sum(dexp(d,rate,log = TRUE))}

EMV2 = mle(NegLogLik2, start = list(rate=1))

## Warning in dexp(d, rate, log = TRUE): NaNs produced

summary(EMV2)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = NegLogLik2, start = list(rate = 1))
## 
## Coefficients:
##       Estimate Std. Error
## rate 0.9562808 0.09562797
## 
## -2 log L: 208.9409

EMV3 = mle2(NegLogLik,start = list(rate=1), data = list(d))

## Warning in dexp(d, rate, log = TRUE): NaNs produced

summary(EMV3)

## Maximum likelihood estimation
## 
## Call:
## mle2(minuslogl = NegLogLik, start = list(rate = 1), data = list(d))
## 
## Coefficients:
##      Estimate Std. Error z value     Pr(z)    
## rate 0.956281   0.095628      10 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## -2 log L: 208.9409

Suponga que cuenta con datos que provienen de una mezcla de dos distribuciones normales, La mayoría de ellos son buenos datos con una varianza pequeña \(\sigma_1^2\) y función de densidad \(f_1\), y una fracción \(p\) de outliers con una varianza mas grande \(σ_2^2\) y función de densidad \(f_2\). Suponga además que ambas medias son iguales a 10. Genere 100 datos aleatorios que cumplan con esta característica y encuentre los EMV para los parámetros \(\mu\), \(\sigma_1^2\), \(\sigma_2^2\) y \(p\), suponiendo que \(\sigma_2^2\) es cuatro veces \(\sigma_1^2\).

set.seed(123)
b = rnorm(n = 100,mean = 10,sd = 0.5)
b

##   [1]  9.719762  9.884911 10.779354 10.035254 10.064644 10.857532 10.230458
##   [8]  9.367469  9.656574  9.777169 10.612041 10.179907 10.200386 10.055341
##  [15]  9.722079 10.893457 10.248925  9.016691 10.350678  9.763604  9.466088
##  [22]  9.891013  9.486998  9.635554  9.687480  9.156653 10.418894 10.076687
##  [29]  9.430932 10.626907 10.213232  9.852464 10.447563 10.439067 10.410791
##  [36] 10.344320 10.276959  9.969044  9.847019  9.809764  9.652647  9.896041
##  [43]  9.367302 11.084478 10.603981  9.438446  9.798558  9.766672 10.389983
##  [50]  9.958315 10.126659  9.985727  9.978565 10.684301  9.887115 10.758235
##  [57]  9.225624 10.292307 10.061927 10.107971 10.189820  9.748838  9.833396
##  [64]  9.490712  9.464104 10.151764 10.224105 10.026502 10.461134 11.025042
##  [71]  9.754484  8.845416 10.502869  9.645400  9.655996 10.512786  9.857613
##  [78]  9.389641 10.090652  9.930554 10.002882 10.192640  9.814670 10.322188
##  [85]  9.889757 10.165891 10.548420 10.217591  9.837034 10.574404 10.496752
##  [92] 10.274198 10.119366  9.686047 10.680326  9.699870 11.093666 10.766305
##  [99]  9.882150  9.486790

set.seed(123)
c = rnorm(n = 100,mean = 10,sd = 2)
c[1]= 15.610303561
c[50]=0.2534628
c[35]=36.23564976524
c[75]=0.0025485
c

##   [1] 15.6103036  9.5396450 13.1174166 10.1410168 10.2585755 13.4301300
##   [7] 10.9218324  7.4698775  8.6262943  9.1086761 12.4481636 10.7196277
##  [13] 10.8015429 10.2213654  8.8883177 13.5738263 10.9957010  6.0667657
##  [19] 11.4027118  9.0544172  7.8643526  9.5640502  7.9479911  8.5422175
##  [25]  8.7499215  6.6266134 11.6755741 10.3067462  7.7237261 12.5076298
##  [31] 10.8529284  9.4098570 11.7902513 11.7562670 36.2356498 11.3772805
##  [37] 11.1078353  9.8761766  9.3880747  9.2390580  8.6105860  9.5841654
##  [43]  7.4692073 14.3379119 12.4159240  7.7537828  9.1942303  9.0666893
##  [49] 11.5599302  0.2534628 10.5066370  9.9429065  9.9142591 12.7372046
##  [55]  9.5484580 13.0329412  6.9024944 11.1692275 10.2477085 10.4318831
##  [61] 10.7592790  8.9953531  9.3335852  7.9628492  7.8564175 10.6070573
##  [67] 10.8964196 10.1060085 11.8445349 14.1001694  9.0179377  5.3816622
##  [73] 12.0114770  8.5815985  0.0025485 12.0511427  9.4304540  7.5585646
##  [79] 10.3626070  9.7222173 10.0115284 10.7705608  9.2586799 11.2887531
##  [85]  9.5590269 10.6635639 12.1936780 10.8703630  9.3481368 12.2976152
##  [91] 11.9870077 11.0967939 10.4774635  8.7441878 12.7213049  8.7994808
##  [97] 14.3746660 13.0652213  9.5285993  7.9471582

boxplot(b, main ="HISTOGRAMA MU=10",xlab = "VALORES DE X", ylab = "FRECUENCIAS", col=5)

boxplot(c, main ="HISTOGRAMA MU=10",xlab = "VALORES DE X", ylab = "FRECUENCIAS", col=5)

NegLogLik = function(mu,sigma){-sum(dnorm(b,mu,sigma,log = TRUE))}

EMV1 = mle(NegLogLik, start = list(mu=10, sigma=0.5))

## Warning in dnorm(b, mu, sigma, log = TRUE): NaNs produced

## Warning in dnorm(b, mu, sigma, log = TRUE): NaNs produced

## Warning in dnorm(b, mu, sigma, log = TRUE): NaNs produced

summary(EMV1)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = NegLogLik, start = list(mu = 10, sigma = 0.5))
## 
## Coefficients:
##         Estimate Std. Error
## mu    10.0452033 0.04541209
## sigma  0.4541209 0.03211057
## 
## -2 log L: 125.909

NegLogLik = function(mu,sigma){-sum(dnorm(c,mu,sigma,log = TRUE))}

EMV1 = mle(NegLogLik, start = list(mu=10, sigma=2))

## Warning in dnorm(c, mu, sigma, log = TRUE): NaNs produced

summary(EMV1)

## Maximum likelihood estimation
## 
## Call:
## mle(minuslogl = NegLogLik, start = list(mu = 10, sigma = 2))
## 
## Coefficients:
##        Estimate Std. Error
## mu    10.312040  0.3510474
## sigma  3.510474  0.2482278
## 
## -2 log L: 534.938

MÁXIMA VEROSIMILITUD

Wilson Sandoval Rodriguez

2022-03-17

Estimación por el Método de Máxima Verosimilitud

Función mle del paquete stats4

Ejercicios