## Ref: RegLiBoot.r ## ############################################################################## ## PROGRAMA PARA ESTIMAR A PRECISAO DOS ESTIMADORES BETA DA REGRESSAO ## ## VIA BOOTSTRAP - Assumindo que o modelo de regressao nao e correto ## ############################################################################## ## Este procedimento seleciona os elementos aleatoriamente, como um ## ## coeficiente de correlacao Bootstrap ## bootstrap=function(n, dados) { ## n e o numero de amostras por bootstrap a partir da amostra inicial e ## ## dados e a amostra inicial ## if(ncol(dados) != 2) stop('dados deve ser uma matrix de 2 colunas') ## se o numero de colunas for diferente de 2 nao se ajusta a regressao ## ## linear simples ## x <- 1:nrow(dados) ## x e um vetor que contem o numero de pares (x,y) na amostra (em ordem ## ## crescente) ## z <-list() ## z e uma lista de objetos a ser construida ## estat <- coef(lm(dados[,2]~dados[,1])) ## estat contem os valores dos coeficientes beta0 e beta1 da amostra ## ## inicial x ## z$estat <-estat ## z$estat e um objeto da lista z que contem os valores dos coeficientes ## ## beta0 e beta1 damostra inicial ¨dados¨. ## data<-matrix(sample(x, size=length(x)*n, replace=T), nrow=n) ## data é uma matriz de n linhas e tantas colunas quantas forem o tamanho ## ## do vetor x, onde cada linha é uma amostra gerada por bootstrap do vetor ## ## x, cujos elementos correspondem as linhas da matriz ¨dados¨ ## bdcoef<-matrix(0, nrow=n,ncol=2) ## bdcoef e uma matriz de n linhas e 2 colunas com zeros como elementos ## ## iniciais. Ou seja, em outras palavras, criamos uma matriz denominada ## for(i in 1:n) { y <- data[i,] bdcoef[i,] <- coef(lm(dados[y,2] ~ dados[y,1])) } ## bdcoef, agora, e uma matriz n linhas e 2 colunas cujos elementos são os ## ## valores dos coeficientes beta p/ cada amostra gerada por bootstrap ## z$distn<-bdcoef ## z$distn e um objeto da lista z que representa o vetor bd ## b0 <- mean(bdcoef[,1]) b1 <- mean(bdcoef[,2]) b <- c(b0,b1) bias<- b-estat z$bias <- bias ## z$bias e um objeto da lista z que contem o vicio, media dos coeficientes ## ## betas calculados p/ cada amostra gerada por bootstrap, em relacao aos ## ## coeficientes calculados sobre a amostra inicial ¨dados¨ ## z$se<-matrix(0, nrow=n,ncol=2) s0=sqrt(var(bdcoef[,1])) s1=sqrt(var(bdcoef[,2])) s <- c(s0,s1) z$se <- s ## z$se e um objeto da lista z que contem o erro padrao dos coeficientes ## ## betas da amostra inicial ¨dados¨ calculada por bootstrap ## z ## retorna a lista construida ## } ## FIM DA FUNCAO BOOTSTRAP ## dados1=read.table('DadRegLiPonto.txt',head=T) ## lendo os dados completos ## x=dados1$PH y=dados1$SB ## extraindo de dados 1 as colunas de interesse: x e y ## dados <- data.frame(x,y) dados ## formando a matriz de dados para regressao linear ## n=200 w <- bootstrap(n, dados) w ############################################################################## ## GRAFICO DA FUNCAO DISTRIBUICAO EMPIRICA DA AMOSTRA INICIAL ## ############################################################################## p=NULL ## cria 1 vetor p que posteriormente contera os valores das funcoes dist. ## ## empiricas de beta0 e beta1 ## par(mfrow=c(2,3)) ## abre uma janela que conterá sete gráficos (3 numa linha e 3 na outra) ## for (i in 1:n) p[i]=i/n plot(sort(w$distn[,1]), p, xlab="Beta0", ylab="F(x) empirica") plot(sort(w$distn[,2]), p, xlab="Beta1", ylab="F(x) empirica") ############################################################################## ## GRAFICO DE DISPERSAO E AJUSTE DO MODELO INICIAL ## ############################################################################## plot(dados, xlab="x", ylab="y") abline(w$estat) ############################################################################## ## GRAFICO DE DISPERSAO E AJUSTES POR BOOTSTRAP ## ############################################################################## plot(dados, xlab="x", ylab="y") for (i in 1:n) abline(w$distn[i,]) ############################################################################## ## HISTOGRAMA DOS VALORES DAS ESTATISTICAS BETA0 E BETA1 ## ############################################################################## hist(w$distn[,1], xlab='Beta0',ylab='Frequencia',main='Distribuicao Bootstrap de Beta0') hist(w$distn[,2], xlab='Beta1',ylab='Frequencia',main='Distribuicao Bootstrap de Beta1') ############################################################################## ## Ref: RegLiBoot.r ## ############################################################################## ## PROGRAMA PARA ESTIMAR A PRECISAO DOS ESTIMADORES BETA DA REGRESSAO ## ## VIA BOOTSTRAP - Assumindo que o modelo de regressao e correto ## ############################################################################## ## Este procedimento ajusta uma regressao, estima os residuos e valores ## ## preditos, e gera amostras bootstrap adicionando aleatoriamente os ## ## residuos amostrados aos valores preditos ## bootstrap=function(n, dados) { ## n e o numero de amostras por bootstrap a partir da amostra inicial e ## ## dados e a amostra inicial ## if(ncol(dados) != 2) stop('dados deve ser uma matrix de 2 colunas') ## se o numero de colunas for diferente de 2 nao se ajusta a regressao ## ## linear simples ## # k <- 1:nrow(dados) ## x e um vetor que contem o numero de pares (x,y) na amostra (em ordem ## ## crescente) ## z <-list() ## z e uma lista de objetos a ser construida ## lm.ini<-lm(dados[,2]~dados[,1]) ## lm.ini e um objeto que contem os coeficiente do modelo, residuos entre ## ## outros estat<-lm.ini$coef ## estat armazena os valores dos coeficientes beta0 e beta1 da amostra ## ## inicial x ## res.ini<-lm.ini$residual ## res.ini armazena os valores dos residuos correspondente a amostra ## ## inicial ## z$estat <-estat ## z$estat e um objeto da lista z que contem os valores dos coeficientes ## ## beta0 e beta1 da amostra inicial ¨dados¨. ## data<-matrix(sample(res.ini, size=length(res.ini)*n, replace=T), nrow=n) ## data e uma matriz de n linhas e tantas colunas quantas forem o tamanho ## ## da amostra inicial, onde cada linha é uma amostra gerada por bootstrap ## ## do vetor que contem os residuos ## ypred=matrix(0, nrow=ncol(data), ncol=n) ## ypred e uma matriz de zeros, com numero de linhas igual ao numero de ## ## elementos da amostra inicial e n colunas ## lmboot=matrix(0, nrow=n, ncol=2) ## lmboot e uma matriz de n linhas e 2 colunas com zeros como elementos ## ## iniciais ## for(i in 1:n) { ypred[,i] = data[i,]+dados[,2] lmboot[i,] <- lm(ypred[,i] ~ dados[,1])$coef } ## lmboot, agora, e uma matriz de n linhas e 2 colunas cujos elementos sao ## ## os valores dos coeficientes beta p/ cada amostra gerada por bootstrap ## z$distn<-lmboot ## z$distn e um objeto da lista z que representa a matriz dos coeficientes ## bias0 <- mean(lmboot[,1])-estat[1] bias1 <- mean(lmboot[,2])-estat[2] ## bias0 contem o valor do vies para beta0 ## ## bias1 contem o valor do vies para beta1 ## z$bias <-c(bias0,bias1) ## z$bias e um objeto da lista z que contem o vicio, media dos coeficientes ## ## betas calculados p/ cada amostra gerada por bootstrap, em relacao aos ## ## coeficientes calculados sobre a amostra inicial ¨dados¨ ## z$se<-matrix(0, nrow=n,ncol=2) ## z$se e uma matriz de n linhas e 2 colunas com zeros como elementos ## ## iniciais s0=sqrt(var(lmboot[,1])) s1=sqrt(var(lmboot[,2])) ## s0 contem o valor do erro padrao para beta0 ## ## s1 contem o valor do erro padrao para beta1 ## z$se <-c(s0,s1) ## z$se e um objeto da lista z que contem o erro padrao dos coeficientes ## ## betas da amostra inicial ¨dados¨ calculada por bootstrap ## z ## retorna a lista construida ## } ## FIM DA FUNCAO BOOTSTRAP ## ############################################################################## dados1=read.table('DadRegLiPonto.txt',head=T) ## lendo os dados completos ## x=dados1$PH y=dados1$SB ## extraindo de dados 1 as colunas de interesse: x e y ## dados <- data.frame(x,y) dados[1:3,] ## formando a matriz de dados para regressao linear ## ## EXECUTANDO A FUNCAO BOOTSTRAP: ## n=200 w <- bootstrap(n, dados) w ############################################################################## ## GRAFICO DA FUNCAO DISTRIBUICAO EMPIRICA DA AMOSTRA INICIAL ## ############################################################################## p=NULL ## cria 1 vetor p que posteriormente contera os valores das funcoes dist. ## ## empiricas de beta0 e beta1 ## par(mfrow=c(2,3)) ## abre uma janela que conterá sete gráficos (4 numa linha e 3 na outra) ## for (i in 1:n) p[i]=i/n plot(sort(w$distn[,1]), p, xlab="Beta0", ylab="F(x) empirica") plot(sort(w$distn[,2]), p, xlab="Beta1", ylab="F(x) empirica") ############################################################################## ## GRAFICO DE DISPERSAO E AJUSTE DO MODELO INICIAL ## ############################################################################## plot(dados, xlab="x", ylab="y") abline(w$estat) ############################################################################## ## GRAFICO DE DISPERSAO E AJUSTES POR BOOTSTRAP ## ############################################################################## plot(dados, xlab="x", ylab="y") for (i in 1:n) abline(w$distn[i,]) ############################################################################## ## HISTOGRAMA DOS VALORES DAS ESTATISTICAS BETA0 E BETA1 ## ############################################################################## hist(w$distn[,1], xlab='Beta0',ylab='Frequencia',main='Distribuicao Bootstrap de Beta0') hist(w$distn[,2], xlab='Beta1',ylab='Frequencia',main='Distribuicao Bootstrap de Beta1')