#
#   amostra_estrat.R - A set of R functions for studying sampling
#                      designs
#
#Copyright (C) 2004 Fernando Henrique F. P. Rosa <feferraz at ime.usp.br>
#
# 
#This program is free software; you can redistribute it and/or
#modify it under the terms of version 2 of the GNU General Public License
#as published by the Free Software Foundation.  A copy of this license should
#be included in the file COPYING.
#
#This program is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#GNU General Public License for more details.
#
#You should have received a copy of the GNU General Public License
#along with this program; if not, write to the Free Software
#Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
#
#

## Esse conjunto de funções em R implementa uma classe 'amostra', para
## estudar diferentes tipos de delineamentos amostrais
## mais detalhes sobre os tipos de delineamentos e exemplos de uso de cada
## função se encontram aqui:
## http://www.linux.ime.usp.br/~feferraz/files/lista/mae315_projeto.pdf


AES <- function(N,nh,variavel,estratos,estatistica=c("media","tau")) {
	Est <- numeric(N)
        Nh <- tapply(as.numeric(estratos),estratos,length)
	for (j in 1:N) {
		ybh <- numeric(max(as.numeric(estratos)))
		for (i in 1:max(as.numeric(estratos))) {
			ybh[i] <- mean(sample(variavel[as.numeric(estratos) == i],nh[i],replace=TRUE))
		}
		Est[j] <- switch(estatistica,media =sum(Nh * ybh)/sum(Nh),tau = sum(Nh * ybh))
	}
	res <- list(resultados=Est,estatistica=estatistica,N=N,n=nh,y=variavel,y.name=deparse(substitute(variavel)),estratos=estratos,estratos.names=deparse(substitute(estratos)),tipo.amostra="AES")
	class(res) <- c("amostra", "list")
        return(res)
}

AASc <- function(N,n,variavel,estatistica=c("media","tau")) {
	estatistica <- match.arg(estatistica)
	Est <- numeric(N)
	for(j in 1:N) {
		yb <- mean(sample(variavel,n,replace=TRUE))
		Est[j] <- switch(estatistica,media = yb,tau = length(variavel)*yb)
	}
        res <- list(resultados=Est,estatistica=estatistica,N=N,n=n,y=variavel,y.name=deparse(substitute(variavel)),tipo.amostra="AASc")
        class(res) <- c("amostra", "list")
        return(res)
}

 

RAAS <- function(N,n,x,y,estatistica=c("razao","tau","media")) {
	if (length(x) != length(y)) { stop("x and y are not the same size.") }
	Est <- numeric(N)
	for (i in 1:N) {
		amostra <- sample(1:length(x),n,replace=TRUE)
		r <- mean(y[amostra])/mean(x[amostra])
		Est[i] <- switch(estatistica,razao = r,tau = r * sum(x),media = r * mean(x))
	}
        res <- list(resultados=Est,estatistica=estatistica,N=N,n=n,y=y,y.name=deparse(substitute(y)),x=x,x.name=deparse(substitute(x)),tipo.amostra="RAAS")
        class(res) <- c("amostra", "list")
        return(res)

}

RAAS.desig <-  function(N,n,x,y,estatistica=c("razao","tau","media")) {
        if (length(x) != length(y)) { stop("x and y are not the same size.") }
        Est <- numeric(N)
        for (i in 1:N) {
                amostra <- sample(1:length(x),n,replace=TRUE,prob=x/sum(x))
                r <- sum(y[amostra]/x[amostra])/n
                Est[i] <- switch(estatistica,razao = r,tau = r * sum(x),media = r * mean(x))
        }
        res <- list(resultados=Est,estatistica=estatistica,N=N,n=n,y=y,y.name=deparse(substitute(y)),x=x,x.name=deparse(substitute(x)),tipo.amostra="RAAS.desig")
        class(res) <- c("amostra", "list")
        return(res)
}

RAAS.estr <-  function(N,nh,x,y,estratos,estatistica=c("tau","media")) { 
        if (length(x) != length(y)) { stop("x and y are not the same size.") }
        Est <- numeric(N) 
	Nh <- tapply(as.numeric(estratos),estratos,length)
        for (i in 1:N) { 
		yrh <-  numeric(max(as.numeric(estratos)))
                for (j in 1:max(as.numeric(estratos))) {
			amostra <- sample(1:Nh[j],nh[j],replace=TRUE) 
			fatores <- as.numeric(estratos)
	                r <- mean(y[fatores == j][amostra])/mean(x[fatores == j][amostra]) 
			yrh[j] <- r * mean(x[fatores == j])
		}
        	Est[i] <- switch(estatistica,tau = sum(Nh*yrh),media = sum(Nh*yrh)/sum(Nh))
        } 
        res <- list(resultados=Est,estatistica=estatistica,estratos=estratos,estratos.name=deparse(substitute(estratos)),N=N,n=nh,y=y,y.name=deparse(substitute(y)),x=x,x.name=deparse(substitute(x)),tipo.amostra="RAAS.estr")
        class(res) <- c("amostra", "list") 
	return(res)
} 

RAAS.estr.desig <-  function(N,nh,x,y,estratos,estatistica=c("tau","media")) {
        if (length(x) != length(y)) { stop("x and y are not the same size.") }
        Est <- numeric(N)
        Nh <- tapply(as.numeric(estratos),estratos,length)
        for (i in 1:N) {
                yrh <-  numeric(max(as.numeric(estratos)))
                for (j in 1:max(as.numeric(estratos))) {
                        fatores <- as.numeric(estratos)
                        amostra <- sample(1:Nh[j],nh[j],replace=TRUE,prob=x[fatores == j]/sum(x[fatores == j]))
                        r <- mean(y[fatores == j][amostra])/mean(x[fatores == j][amostra])
                        yrh[j] <- r * mean(x[fatores == j])
                }
                Est[i] <- switch(estatistica,tau = sum(Nh*yrh),media = sum(Nh*yrh)/sum(Nh))
        }
        res <- list(resultados=Est,estatistica=estatistica,estratos=estratos,estratos.name=deparse(substitute(estratos)),N=N,n=nh,y=y,y.name=deparse(substitute(y)),x=x,x.name=deparse(substitute(x)),tipo.amostra="RAAS.estr.desig")
        class(res) <- c("amostra", "list")
        return(res)
}


print.amostra <- function(res) {
	cat("    Análise do Procedimento Amostral\n\n")
	cat("  Procedimento amostral: ",res$tipo.amostra,"\n",sep='')
	cat("  Variavel: ",res$y.name,"\n")
	cat("  N = ",length(res$y),", n = ",sum(res$n),"\n",sep='')
	if (!is.null(res$estratos)) {
        cat("  estratos: ",res$estratos.name," (",length(res$n)," níveis)\n",sep='')
	cat("  nh:",res$n,"\n")
	}
	if(!is.null(res$x)) {
	cat("  Variável auxiliar X: ",res$x.name,"\n")
	}
	cat("  Numero de amostras simuladas: ",res$N,"\n")
	cat("  Estatistica utilizada: ",res$estatistica,"\n")
	cat("\tvalor populacional: ",switch(res$estatistica,tau = sum(res$y),media = mean(res$y), razao = sum(res$y)/sum(res$x)),"\n",sep='') 
	cat("\tmedia: ",format(mean(res$resultados),nsmall=0),"\n")
	cat("\tvariância: ",var(res$resultados),"\n")
	return(invisible(res))
}

plot.amostra <- function(res,qqplot=FALSE,...) {
	if (qqplot == TRUE) { par(mfrow=c(1,2)) }

	titulo <- paste(res$tipo.amostra,"\n",res$N,"amostras de tamanho",sum(res$n))
	xlab <- paste(res$y.name,"\nEstatistica: ",res$estatistica)
	parametro <- switch(res$estatistica,tau = sum(res$y),media = mean(res$y), razao = sum(res$y)/sum(res$x))

	hist(res$resultados,main=titulo,xlab=xlab,prob=T,...)
	abline(v=parametro,lty=2)

	if (qqplot == TRUE) { qqnorm(res$resultados) }
}

compara <- function() {
# graficos da comparacao

aes <- AES(100,nh.otima,celulares,regiao,estatistica="tau")
aas <- AASc(1000,10,celulares,estatistica="tau")
raas <- RAAS(100,10,populacao,celulares,estatistica="tau")#
raas.estr <- RAAS.estr(100,nhr.otimo,populacao,celulares,regiao,estatistica="tau")
#pdf('expansao.pdf')
par(mfrow=c(2,2))
minimo <- min(aes,aas,raas,raas.estr)
maximo <- max(aes,aas,raas,raas.estr)
xlimite <- c(minimo,maximo)
hist(aas,xlim=xlimite,prob=T)
hist(aes,xlim=xlimite,prob=T)
hist(raas,xlim=xlimite,prob=T)
hist(raas.estr,xlim=xlimite,prob=T)


hist(aas,xlim=c(min(aas,aes),max(aas,aes)),prob=T,main="Estimador Expansao AASc (n=10)",xlab="T") 
abline(v=sum(celulares),lty=2)
hist(aes,xlim=c(min(aas,aes),max(aas,aes)),prob=T,main="Estimador Expansao AES (n=10)",xlab="Tes")
abline(v=sum(celulares),lty=2)
dev.off()


}