# Alejandro Pozo Pazos 
# Boletín 2 - Ejercicio 2
# Distribución logística. 
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#

# Gráfica con las funciones de densidad de Cauchy,
# Normal y Logística de parámetros 0 y 1.
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#

# Gráfica con las funciones de densidad de Cauchy, Normal y Logística
# de parámetros 0 y 1.
curve(dcauchy(x,0,1),axes=FALSE,-4,4,ylab="",xlab="",col="green",lty=3,lwd=3,ylim=c(0,0.5))
curve(dlogis(x,0,1), -4,4,add=TRUE,col="red",lty=1,lwd=3)
curve(dnorm(x,0,1),-4,4,add=TRUE,col="blue",lty=5,lwd=3)
axis(1,pos=0)
axis(2,pos=-4)
legend(1.5, 0.5, c("Logística", "Normal", "Cauchy"),col=c("red","blue","green"),lty=c(3,1,5))
title(sub="Funciones de densidad de Cauchy, normal y logística de parámetros (0,1)",col.sub="gray10",cex.sub=1,font.sub=1,line=2)

# Generador de valores para la distribución logística.
# Método de inversión.
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#

log_inv<-function(n,alpha=0,beta=1){
	t<-proc.time()
	x<-runif(n,0,1)
	x<-alpha+beta*log(x/(1-x))
	t<-proc.time()-t
	return(list("datos"=x,"tiempo"=t))
}

# Contraste de Kolmogorov-Smirnov.
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#

TEST_KS<-function(x,...){
	result<-ks.test(x, ...)
	sal<-c(result$p.value,result$statistic)
	return(sal)
}

# Programa para evaluar la bondad de ajuste del generador
# creado con el método de inversión.
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#

# Establecemos los parámetros
alpha=0
beta=1

# Generamos las secuencias de valores
K=100
N=10000
x<-log_inv(N*K,alpha,beta)
cat("Tiempo de generación:",x$tiempo[[3]],"segundos.\n")
x<-x$datos
# Convertimos x en una matriz de 100 columnas
dim(x)<-c(N,K)

# Por columnas, realizamos el contraste de KS
res_test<-apply(x,2,TEST_KS,"plogis",alpha,beta)

# Creamos un vector con los p-valores 
# y otro con los estadísticos del contraste
(vector_pvalores<-res_test[1,])
(vector_estadistico<-res_test[2,])

# Nos quedamos con los que se rechazan:
rechazados<-vector_pvalores[vector_pvalores<0.05]

# Devolvemos los resultados
(pvalor_estimado<-length(rechazados)/K)
summary(vector_estadistico)

# Generador de valores para la distribución logística.
# Método de aceptación rechazo, a partir de Cauchy o de Laplace.
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#

# Función de densidad de Laplace
dlaplace<-function(x,alpha=1){
	dens<-(alpha/2)*exp(-alpha*abs(x))
	return(dens)
}

# Gráficas de los cocientes Logística/Cauchy y Logística/Laplace
op<-par(mfrow=c(2,1))
curve(dlogis(x,0,1)/dcauchy(x,0,1),-10,10,
ylim=c(0,2.5),lwd=2,ylab="",xlab="",col="green",
main="Logístiva vs. Cauchy")
abline(h=2, lty=2,lwd=2)
curve(dlogis(x,0,1)/dlaplace(x,1),-10,10,
ylim=c(0,2.5),lwd=2,ylab="",xlab="",col="red",
main="Logístiva vs. Laplace")
abline(h=2, lty=2,lwd=2)
par(op)

# Generar Cauchy por el método de inversión
cauchy<-function(n,sigma){
	x<-runif(n,0,1)
	x<-sigma*tan(pi*(x-0.5))
	return(x)
}

# Generar Laplace por el método de convolución
laplace<-function(n,alpha=1){
	x<-runif(n,0,1)
	x<-(-1/alpha)*log(x)
	u<-runif(n)
	x[u>0.5]<-(-1.0)*x[u>0.5]
	return(x)
}

# Generar logística
log_ar<-function(n,op){
	x<-numeric(n)
	int<-0
	ngen<-0
	t<-proc.time()
	if (op=="cauchy"){
		c<-2
		while(ngen<n){
			int<-int+1
			v<-cauchy(1,1)
			u<-runif(1)*c*dcauchy(v,0,1)
			if (u<dlogis(v,0,1)) {
				ngen<-ngen+1
				x[ngen]<-v
			}
		}
	} else if (op=="laplace") {
		c<-2
		while(ngen<n){
			int<-int+1
			v<-laplace(1,1)
			u<-runif(1)*c*dlaplace(v,1)
			if (u<dlogis(v)) {
				ngen<-ngen+1
				x[ngen]<-v
			}
		}
	}
	t=proc.time()-t
	return(list("datos"=x,"intentos"=int,"tiempo"=t))
}

# Programa para evaluar la bondad de ajuste del generador
# creado con el método de aceptación-rechazo.
#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#-#

# Generamos las secuencias de valores
K=100
N=1000
x1<-log_ar(N*K,"cauchy")
x2<-log_ar(N*K,"laplace")
cat("Tiempo de generación con Cauchy:",x1$tiempo[[3]],"segundos.\n",
"Tiempo de generación con Laplace:",x2$tiempo[[3]],"segundos.\n")
x1<-x1$datos
x2<-x2$datos

# Convertimos x en una matriz de 100 columnas
dim(x1)<-c(N,K)
dim(x2)<-c(N,K) 

# Por columnas, realizamos el contraste de KS
res_test1<-apply(x1,2,TEST_KS,"plogis",0,1)
res_test2<-apply(x2,2,TEST_KS,"plogis",0,1)

# Creamos un vector con los p-valores 
# y otro con los estadísticos del contraste
(vector_pvalores1<-res_test1[1,])
(vector_pvalores2<-res_test2[1,])
(vector_estadistico1<-res_test1[2,])
(vector_estadistico2<-res_test2[2,])

# Nos quedamos con los que se rechazan:
rechazados1<-vector_pvalores1[vector_pvalores1<0.05]
rechazados2<-vector_pvalores2[vector_pvalores2<0.05]

# Devolvemos los resultados
(pvalor_estimado1<-length(rechazados1)/K)
(pvalor_estimado2<-length(rechazados2)/K)
summary(vector_estadistico1)
summary(vector_estadistico2)