Rut:19.474.068-9

Link:http://rpubs.com/JoseTapiaJara/391377

Optimizacion de máxima vero-similitud utilizando la función “optim”

Primera parte: importar datos

Se dispone de un archivo llamado “empleados”, el cual contiene los datos para comparar este taller y será importado por medio del código:
library(readxl)
empleados <- read_excel("C:/Users/cotel/Documents/R/empleados.xls")
empleados$genero  <- factor(empleados$genero, labels=c("Femenino", "Masculino"))
empleados$catlab  <- factor(empleados$catlab, labels=c("Administrativo", "Seguridad", "Directivo"))
empleados$minoria <- factor(empleados$minoria, labels=c("No", "Si"))
#Visualiza los primeros datos
head(empleados)
## # A tibble: 6 x 9
##      id genero     educ catlab     salario salini tiempemp expprev minoria
##   <dbl> <fct>     <dbl> <fct>        <dbl>  <dbl>    <dbl>   <dbl> <fct>  
## 1     1 Masculino    15 Directivo    57000  27000       98     144 No     
## 2     2 Masculino    16 Administr~   40200  18750       98      36 No     
## 3     3 Femenino     12 Administr~   21450  12000       98     381 No     
## 4     4 Femenino      8 Administr~   21900  13200       98     190 No     
## 5     5 Masculino    15 Administr~   45000  21000       98     138 No     
## 6     6 Masculino    15 Administr~   32100  13500       98      67 No
# Estructura del dataframe
str(empleados)    
## Classes 'tbl_df', 'tbl' and 'data.frame':    474 obs. of  9 variables:
##  $ id      : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ genero  : Factor w/ 2 levels "Femenino","Masculino": 2 2 1 1 2 2 2 1 1 1 ...
##  $ educ    : num  15 16 12 8 15 15 15 12 15 12 ...
##  $ catlab  : Factor w/ 3 levels "Administrativo",..: 3 1 1 1 1 1 1 1 1 1 ...
##  $ salario : num  57000 40200 21450 21900 45000 ...
##  $ salini  : num  27000 18750 12000 13200 21000 ...
##  $ tiempemp: num  98 98 98 98 98 98 98 98 98 98 ...
##  $ expprev : num  144 36 381 190 138 67 114 0 115 244 ...
##  $ minoria : Factor w/ 2 levels "No","Si": 1 1 1 1 1 1 1 1 1 1 ...
# Nombre de las variables
names(empleados)   
## [1] "id"       "genero"   "educ"     "catlab"   "salario"  "salini"  
## [7] "tiempemp" "expprev"  "minoria"
# Separa los nombre de cada variable del dataframe
attach(empleados)

# Resumen de de una variable 
mean(salario)
## [1] 34419.57
median(salario)
## [1] 28875
var(salario) # Variance
## [1] 291578214
sd(salario) # Standard deviation
## [1] 17075.66
max(salario) # Max value
## [1] 135000
min(salario) # Min value
## [1] 15750
range(salario) # Range
## [1]  15750 135000
quantile(salario) # Quantiles 25%
##       0%      25%      50%      75%     100% 
##  15750.0  24000.0  28875.0  36937.5 135000.0
quantile(salario, c(.3,.6,.9)) # Customized quantiles
##     30%     60%     90% 
## 24885.0 30750.0 59392.5
Este fragmento de código además de importar dichos datos, nos servirá para explorar sus propiedades(como por ejemplo la media,máximo,mínimo,etc) y además de graficarlos por medio de los siguientes comandos:
summary(salario)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   15750   24000   28875   34420   36938  135000
hist(salario)

boxplot(salario)

stem(salario)
## 
##   The decimal point is 4 digit(s) to the right of the |
## 
##    1 | 666666777777777778888999
##    2 | 00000000000000111111111111111111122222222222222222222222233333333333+148
##    3 | 00000000000000000001111111111111111111111111122222222222223333333333+36
##    4 | 0000000001112222334445555666778899
##    5 | 0111123344555556677778999
##    6 | 0001122355566777888999
##    7 | 00134455889
##    8 | 01346
##    9 | 1127
##   10 | 044
##   11 | 1
##   12 | 
##   13 | 5
Y podemos ordenarlo de la siguiente manera:
resumengrupo <- tapply(salario,genero, summary) 
resumengrupo
## $Femenino
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   15750   21563   24300   26032   28500   58125 
## 
## $Masculino
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   19650   28050   32850   41442   50413  135000

Segunda Parte: Identificar la variable Salario

Esto será por medio de los parametros anteriores, gracias a esta variable será posible buscar algun modelo de distribución matemática(alguna función con diversos parametros) que se ajuste a estos datos, gracias a esto podremos optimizar dicha función con el método “optim” para identificar el valor de dichos parametros y graficarla para tener una visión mucho mas tangible de esto.

Tercera Parte: Modelar

por medio de la instrucción
hist(salario)
hist(exp())
es posible modelar los datos y recibir un histograma, luego de esto compararlos con el histograma de la funcion exponencial nos da paso a establecer un parecido entre esta función y el modelo exponencial, por lo tanto:
fde<-function(x,teta){
  return(teta*exp(-teta*x))
}
Definiremos fde(x, teta) como nuestro modelo exponencial

Cuarta Parte: Definir la funcion Log-Verosimilitud

Esta es la manera en la cual optimizaremos la función mediante el uso de la formula entregada por el profesor. Para optimizar una función es equivalente a optimizar su logaritmo, entonces:
logf<-function(p,x){
  n=length(x)
  f=c()
  for( i in 1:n) {
    f=c(f,-(log(p)-p*x[i]))
  }
  sumaf=sum(f)
  return(sumaf)
}

Quinta Parte: Optimizar la función de Log-F

Esto lo haremos utilizando el método optim, asignandole el resultado a una variable para obtener el valor de tetha:
resultado = optim(par=1/mean(salario),fn=logf,x=salario)
## Warning in optim(par = 1/mean(salario), fn = logf, x = salario): one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
print(resultado)
## $par
## [1] 2.905324e-05
## 
## $value
## [1] 5425.584
## 
## $counts
## function gradient 
##       18       NA 
## 
## $convergence
## [1] 0
## 
## $message
## NULL
tetha= resultado$par[1]
print(tetha)
## [1] 2.905324e-05

Sexta Parte:

Luego de estimado el parametro tetha, damos paso a graficar el histograma de Salario para, de esta forma, compararlo con la curva entregada del modelo exponencial.
hist(salario,freq = F)
curve(fde(x,tetha), add = T)

Al comparar ambos modelos queda en evidencia que mientras los salarios más suben(a la derecha del gráfico), esto más se acerca al modelo exponencial con el valor encontrado de tetha, pero cuando nos acercamos al extremo izquierdo del gráfico nos encontramos con que el primer fragmento del histograma es menor que la curva dada, y el siguiente es mayor que dicha función modelada. Pese a que se asemeja mucho, esos son los detalles que se observan al comparar ambos gráficos.