# Análisis de Datos.



library(mea, lib.loc="../")


# Estudio de Simetría.


# Gráfico de Simetría.

simetria.plot<-function(x,...){
a<-sort(x)
b<-(-sort(-x))
plot((b-a)[(b-a)>=0],(b+a)[(b-a)>=0]/2,
     main="Gráfico de Simetría", xlab="x",ylab="y",...)}



# En el archivo faithful tenemos un data.frame con dos variables
# (dos columnas): 
#   1. eruptions, que corresponde a las duraciones de cada erupción en minutos,
#   2. waiting, que corresponde a los tiempos de espera, en minutos, entre dos erupciones,
# en ambos casos para el geyser "Old Faithful" del Parque Nacional de Yellowstone (Wyoming, 
# USA).

data(faithful)
attach(faithful)

stem(eruptions)
hist(eruptions,br=20)


x<-eruptions[eruptions > 3]

stem(x)
hist(x)

simetria.plot(x,col=2)



# Test de Hipótesis para el Coef. de Asimetría

simetria.test<-function(x,alfa=0.05){
 n<-length(x)
 if(n < 50) stop("la muestra es de tamaño insuficiente: n < 50")
 gamma1<-mean((x-mean(x))^3)/(sd(x)^3)
 IC<-c(gamma1-qnorm(1-alfa/2)*sqrt(6/n),gamma1+qnorm(1-alfa/2)*sqrt(6/n))
 EC<-abs(gamma1/sqrt(6/n))
 p<-(1-pnorm(EC))*2
 cat("\n Estudio de simetría. H0: los datos proceden de una distribución simétrica \n\n")
 cat("Coef. de asimetría: ", gamma1, "\n")
 cat("IC. para el coef. de asimetría, al nivel ",1-alfa,": (",
      IC[1], "," ,IC[2],") \n")
 cat("p-valor: ", p ,"\n\n")
}


simetria.test(x)




simetria.plot(x^2,col=2)
simetria.test(x^2)


simetria.plot(x^3,col=2)
simetria.test(x^3)




y<-eruptions[eruptions <= 3]

stem(y,2)


simetria.plot(y,col=2)
simetria.test(y)


simetria.plot(sqrt(y),col=2)
simetria.test(sqrt(y))

simetria.plot(log(y),col=2)
simetria.test(log(y))

simetria.plot(-1/y,col=2)
simetria.test(-1/y)



# Transformaciones de Box-Cox

boxcox<-function(x,alfa){ # Consideramos x un vector positivo
  if(alfa==0) log(x)
  else (x^alfa-1)/alfa
}


simetria.plot(boxcox(y,-1),col=2)
simetria.test(boxcox(y,-1))

simetria.plot(boxcox(y,-2),col=2)
simetria.test(boxcox(y,-2))

simetria.plot(boxcox(y,-3),col=2)
simetria.test(boxcox(y,-3))




# Estudio de Normalidad.


# Q-Q Plots.

qqnorm(x);qqline(x,col=2)


par(mfrow=c(2,2))

qqnorm(x);qqline(x,col=2)
for(i in 1:3) {qqnorm(z<-rnorm(length(x),mean(x),sd(x)),main="Q-Q.Simulación")
               qqline(z,col=2)}

par(mfrow=c(1,1))



# Estimación de la Densidad.


hist(x,prob=T,col=5)
lines(density(x),col=2)
lines(z<-seq(min(x),max(x),length=1000),dnorm(z,mean(x),sd(x)),col=4)



# Función de Distribución Empírica.


ecdf.cont(x)
lines(z<-seq(min(x),max(x),length=1000),pnorm(z,mean(x),sd(x)),col=4)



par(mfrow=c(2,2))

ecdf.cont(x)
lines(z<-seq(min(x),max(x),length=1000),pnorm(z,mean(x),sd(x)),col=4)

for(i in 1:3) {
z<-rnorm(length(x),mean(x),sd(x))
ecdf.cont(z,main="ecdf.simulación")
lines(zz<-seq(min(z),max(z),length=1000),pnorm(zz,mean(x),sd(x)),col=4)
}

par(mfrow=c(1,1))



# Test de Hipótesis para Normalidad.

library(ctest)

ks.test(x,"pnorm",mean(x),sd(x))

shapiro.test(x)



ks.test(x^3,"pnorm",mean(x^3),sd(x^3))

shapiro.test(x^3)



# Estudio de Normalidad para la variable y<-eruptions[eruptions <= 3]

qqnorm(y);qqline(y,col=2)

par(mfrow=c(2,2))

qqnorm(y);qqline(y,col=2)
for(i in 1:3) {qqnorm(z<-rnorm(length(y),mean(y),sd(y)),main="Q-Q.simulación")
               qqline(z,col=2)}

par(mfrow=c(1,1))


ks.test(y,"pnorm",mean(y),sd(y))
shapiro.test(y)


z<-boxcox(y,-3)
ks.test(z,"pnorm",mean(z),sd(z))
shapiro.test(z)




# Estudio de Homocedasticidad.


# En el archivo iris tenemos las medidas en centímetros de las variables:
# 1. Longitud y anchura de los sepalos (Sepal.Length, Sepal.Width, respectivamente).
# 2. Longitud y anchura de los petalos (Petal.Length, Petal.Width, respectivamente).
# de 50 flores de cada una de 3 Especies de lirios (Setosa, Versicolor y Virgínica).
# (Datos publicados por E. Anderson(1935) y R.A. Fisher(1936))


data(iris)
attach(iris)


# Estudiemos si las anchuras de los sepalos se distribuyen con igual
# varianza entre las 3 especies de lirios consideradas.

boxplot(Sepal.Width~Species,col=5)

shapiro.test(Sepal.Width[Species=="setosa"])
shapiro.test(Sepal.Width[Species=="versicolor"])
shapiro.test(Sepal.Width[Species=="virginica"])


bartlett.test(Sepal.Width,Species)

fligner.test(Sepal.Width,Species)


# Estudiemos si las longitudes de los sepalos se distribuyen con igual
# varianza entre las 3 especies de lirios consideradas.

boxplot(Sepal.Length~Species,col=5)

shapiro.test(Sepal.Length[Species=="setosa"])
shapiro.test(Sepal.Length[Species=="versicolor"])
shapiro.test(Sepal.Length[Species=="virginica"])


bartlett.test(Sepal.Length,Species)

fligner.test(Sepal.Length,Species)



boxplot(boxcox(Sepal.Length,-0.8)~Species,col=5)

shapiro.test(boxcox(Sepal.Length[Species=="setosa"],-0.8))
shapiro.test(boxcox(Sepal.Length[Species=="versicolor"],-0.8))
shapiro.test(boxcox(Sepal.Length[Species=="virginica"],-0.8))

bartlett.test(boxcox(Sepal.Length,-0.8),Species)

fligner.test(boxcox(Sepal.Length,-0.8),Species)