Estimacion Puntual
Sebastian Chaparro
21/3/2022
Solucion Taller
Punto 1.
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.
Primero generamos los 50 datos con distribucion Poisson con \(\lambda=5\)
set.seed(123)
DatosPoisson <- rpois(50,5)
GraficaP <- ggplot(data.frame(DatosPoisson),aes(x=DatosPoisson))+geom_bar(fill="orange")+ggtitle("Visualizacion de datos", "Poisson")+ylab("Frecuencia")+xlab("Datos")+theme_cowplot(10)
Graficad <- ggplot(data.frame(DatosPoisson),aes(x=DatosPoisson))+geom_density()+ggtitle("Visualizacion de datos", "Poisson")+ylab("Frecuencia")+xlab("Datos")+theme_cowplot(10)
GraficaP/Graficad
Creamos la funcion de lg-maximaverosimilitud para los datos con dis Poisson
maxveroPoisson <- function(lambda){
-sum(dpois(DatosPoisson,lambda,log = T))
}Método mle
summary(mle(maxveroPoisson, start = 2))## Maximum likelihood estimation
##
## Call:
## mle(minuslogl = maxveroPoisson, start = 2)
##
## Coefficients:
## Estimate Std. Error
## [1,] 5.24 0.3237283
##
## -2 log L: 219.0201Método mle2
summary(mle2(maxveroPoisson,start = list(lambda=5),data = list(DatosPoisson) ))## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = maxveroPoisson, start = list(lambda = 5), data = list(DatosPoisson))
##
## 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.0201Método mle2 sin funcion logmaxverosimilitud
summary(mle2(DatosPoisson~dpois(lambda),start = list(lambda=5) , data= data.frame(DatosPoisson)))## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = DatosPoisson ~ dpois(lambda), start = list(lambda = 5),
## data = data.frame(DatosPoisson))
##
## 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.0201Punto 2.
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.
set.seed(123)
xwibull <- rweibull(50,shape = 3,scale = 200)
Grafica2d <- ggplot(data.frame(xwibull),aes(x=xwibull))+geom_histogram(fill="orange",color="black")+ggtitle("Visualizacion de datos", "Weibull")+ylab("Frecuencia")+xlab("Datos")+theme_cowplot(10)
Grafica1d <- ggplot(data.frame(xwibull),aes(x=xwibull))+geom_density()+ggtitle("Visualizacion de datos", "Weibull")+ylab("Frecuencia")+xlab("Datos")+theme_cowplot(10)
Grafica2d/Grafica1d## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Dis Poisson
Creamos la funcion de lg-maximaverosimilitud para los datos con dis Poisson
maxverowibull<- function(shape,scale){
-sum(dweibull(xwibull,shape,scale,log = T))
}Método mle.
summary(mle(maxverowibull, start = list(shape=1,scale=100)))## Maximum likelihood estimation
##
## Call:
## mle(minuslogl = maxverowibull, start = list(shape = 1, scale = 100))
##
## Coefficients:
## Estimate Std. Error
## shape 2.984199 0.3376886
## scale 193.707655 9.6500462
##
## -2 log L: 556.8388método mle2.
summary(mle2(maxverowibull,start = list(shape=2,scale=99),data = list(xwibull) ))## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = maxverowibull, start = list(shape = 2, scale = 99),
## data = list(xwibull))
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## shape 2.98380 0.33765 8.8369 < 2.2e-16 ***
## scale 193.70524 9.65110 20.0708 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 556.8388método mle2 sin la funcion creada
summary(mle2(xwibull~dweibull(shape,scale),start = list(shape=1,scale=200) , data= data.frame(xwibull)))## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = xwibull ~ dweibull(shape, scale), start = list(shape = 1,
## scale = 200), data = data.frame(xwibull))
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## shape 2.98633 0.33783 8.8398 < 2.2e-16 ***
## scale 193.88857 9.65689 20.0778 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 556.8392Dis Gamma
set.seed(123)
xgamma <- rgamma(n = 50,shape = 3,scale = 3)
graficoGamma <- ggplot(data.frame(xgamma),aes(x=xgamma))+geom_density()+ggtitle("Visualizacion de datos", "Gamma")+ylab("Frecuencia")+xlab("Datos")+theme_cowplot(10)
graficoGamma2 <- ggplot(data.frame(xgamma),aes(x=xgamma))+geom_histogram(fill="orange",color="black")+ggtitle("Visualizacion de datos", "Gamma")+ylab("Frecuencia")+xlab("Datos")+theme_cowplot(10)
graficoGamma2/graficoGamma## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
maxverogamma<- function(shape,scale){
-sum(dgamma(xgamma,shape,scale,log = T))
}Método mle.
summary(mle(maxverogamma, start = list(shape=1,scale=1)))## Maximum likelihood estimation
##
## Call:
## mle(minuslogl = maxverogamma, start = list(shape = 1, scale = 1))
##
## Coefficients:
## Estimate Std. Error
## shape 2.8679241 0.54352580
## scale 0.3859042 0.07992057
##
## -2 log L: 277.1207Método mle2
summary(mle2(maxverogamma,start = list(shape=2,scale=3),data = list(xgamma) ))## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = maxverogamma, start = list(shape = 2, scale = 3),
## data = list(xgamma))
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## shape 2.867949 0.543544 5.2764 1.318e-07 ***
## scale 0.385909 0.079924 4.8285 1.376e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 277.1207Metodo mle2 sin funcion.
summary(mle2(xgamma~dgamma(shape,scale),start = list(shape=1,scale=30) , data= data.frame(xgamma)))## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = xgamma ~ dgamma(shape, scale), start = list(shape = 1,
## scale = 30), data = data.frame(xgamma))
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## shape 2.868241 0.543601 5.2764 1.318e-07 ***
## scale 0.385945 0.079931 4.8285 1.376e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 277.1207Dis Chi-Cuadrado
set.seed(123)
xchi <- rchisq(50,29)
graficaxchi <- ggplot(data.frame(xchi),aes(x=xchi))+geom_histogram(fill="blue",alpha=0.7,color="black")+ggtitle("Visualizacion de datos", "Chi-cuadrado")+ylab("Frecuencia")+xlab("Datos")+theme_cowplot(10)
graficaxchi2 <- ggplot(data.frame(xchi),aes(x=xchi))+geom_density(color="blue")+ggtitle("Visualizacion de datos", "Chi-cuadrado")+ylab("Frecuencia")+xlab("Datos")+theme_cowplot(10)
graficaxchi/graficaxchi2## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Método mle
maxverochi <- function(df){
-sum(dchisq(xchi,df,log = TRUE))
}
summary(mle(maxverochi,start = 10))## Maximum likelihood estimation
##
## Call:
## mle(minuslogl = maxverochi, start = 10)
##
## Coefficients:
## Estimate Std. Error
## [1,] 27.48469 1.02951
##
## -2 log L: 331.1864Método mle2
summary(mle2(maxverochi,start = list(df=2),data = list(xchi) ))## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = maxverochi, start = list(df = 2), data = list(xchi))
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## df 27.4847 1.0295 26.697 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 331.1864Método mle2 sin funcion creada
summary(mle2(xchi~dchisq(df),start = list(df=40) , data= data.frame(xchi)))## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = xchi ~ dchisq(df), start = list(df = 40), data = data.frame(xchi))
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## df 27.4847 1.0295 26.697 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 331.1864EXPONENCIAL
set.seed(123)
xexpo <- rexp(50,2)
plotExpo <- ggplot(data.frame(xexpo),aes(x=xexpo))+geom_density(color="red")+ggtitle("Visualizacion de datos", "Exponencial")+ylab("Frecuencia")+xlab("Datos")+theme_cowplot(10)
plotExpo2 <- ggplot(data.frame(xexpo),aes(x=xexpo))+geom_histogram(fill="red",color="black")+ggtitle("Visualizacion de datos", "Exponencial")+ylab("Frecuencia")+xlab("Datos")+theme_cowplot(10)
plotExpo/plotExpo2 ## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Funcion Maxima verosimilitud.
maxveroExpo <- function(rate){-sum(dexp(xexpo,rate,log = T))}Método mle
summary(mle(maxveroExpo,start = 10))## Warning in dexp(xexpo, rate, log = T): NaNs produced
## Warning in dexp(xexpo, rate, log = T): NaNs produced## Maximum likelihood estimation
##
## Call:
## mle(minuslogl = maxveroExpo, start = 10)
##
## Coefficients:
## Estimate Std. Error
## [1,] 1.769271 0.2502127
##
## -2 log L: 42.93986Método mle2
summary(mle2(maxveroExpo,start = list(rate=2),data = list(xexpo) ))## Warning in dexp(xexpo, rate, log = T): NaNs produced## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = maxveroExpo, start = list(rate = 2), data = list(xexpo))
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## rate 1.76933 0.25022 7.0711 1.537e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 42.93986Método mle sin funcion lg max verosimilitud
summary(mle2(xexpo~dexp(rate),start = list(rate=3) , data= data.frame(xexpo)))## Warning in dexp(x = c(0.421728630529201, 0.288305135443807, 0.664527433903372, :
## NaNs produced## Warning in dexp(x = c(0.421728630529201, 0.288305135443807, 0.664527433903372, :
## NaNs produced## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = xexpo ~ dexp(rate), start = list(rate = 3),
## data = data.frame(xexpo))
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## rate 1.76933 0.25022 7.0711 1.537e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 42.93986Punto 3
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)
DatosPunto3 <- rnorm(100,10,0.5)
DatosPunto3[1]=11.5
DatosPunto3[50]=13
DatosPunto3[100]=5
DatosPunto3[3]=1
DatosPunto3[25]=3.5
DatosPunto31 <- data.frame(DatosPunto3)
ggplot(DatosPunto31,aes(x="",y=DatosPunto3))+geom_boxplot(fill="orange")+
ggtitle("Grafico","BoxPlot") + xlab("Variable") + ylab("Valores")
Funcion Maxima verosimilitud.
MaxVero3 = function(mu,sd){
-sum(dnorm(DatosPunto3,mu,sd,log=T))
}Método mle
summary(mle(MaxVero3,start = list(mu=8,sd=4)))## Maximum likelihood estimation
##
## Call:
## mle(minuslogl = MaxVero3, start = list(mu = 8, sd = 4))
##
## Coefficients:
## Estimate Std. Error
## mu 9.888886 0.13343297
## sd 1.334330 0.09435116
##
## -2 log L: 341.4735Metodo mle2
summary(mle2(MaxVero3,start = list(mu=10,sd=20),data = list(DatosPunto3)))## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = MaxVero3, start = list(mu = 10, sd = 20), data = list(DatosPunto3))
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## mu 9.888885 0.133433 74.111 < 2.2e-16 ***
## sd 1.334329 0.094351 14.142 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 341.4735Método mle sin funcion lg max verosimilitud
summary(mle2(DatosPunto3~dnorm(mu,sd),start = list(mu=19,sd=9) , data= data.frame(DatosPunto3)))## Maximum likelihood estimation
##
## Call:
## mle2(minuslogl = DatosPunto3 ~ dnorm(mu, sd), start = list(mu = 19,
## sd = 9), data = data.frame(DatosPunto3))
##
## Coefficients:
## Estimate Std. Error z value Pr(z)
## mu 9.888892 0.133434 74.111 < 2.2e-16 ***
## sd 1.334336 0.094352 14.142 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## -2 log L: 341.4735