####### Códigos Aula 1 - Inferência 2 ######### ### Exemplo de Binomial ## Vamos simular um p qualquer que será como o parâmetro populacional da Binomial(1,p) set.seed(5000) p <- runif(1,0,1) ## Dado este p populacional vamos retirar uma amostra desta população amostra10 <- rbinom(10,size = 1, prob=p) amostra30 <- rbinom(30,size = 1, prob=p) amostra50 <- rbinom(50,size = 1, prob=p) amostra100 <- rbinom(100,size = 1, prob=p) amostra1000 <- rbinom(1000,size = 1, prob=p) ## Problema geral - Qual é o p ? conjunta <- function(amostra,p){ saida = prod(p^amostra * (1-p)^(1-amostra)) return(saida) } log.conjunta <- function(amostra,p){ saida <- sum(amostra*log(p)) + sum((1-amostra)*log(1-p)) return(saida) } log.conjunta(amostra10,p=0.8) conjunta(amostra10,p=0.8) ## Para testar log(conjunta(amostra10,p=0.8)) exp(log.conjunta(amostra10,p=0.8)) ## Podemos interpretar esse valor como sendo a plausibilidade deste dado p ter gerado aquela particular configuração amostral. ## Sendo assim, é razoável pensar em usar como estimativa de p, o valor mais provável de ter gerado esta particular amostra. ## Mas a distribuição conjunta é uma função da amostra para um dado p parâmetro, sendo assim precisamos olhar a conjunta não #mais em função da amostra que já foi observada e neste momento é uma quantidade fixa, mais sim em função de p que não sabemos # quanto é, assim temos verossimilhanca <- function(p,amostra){ saida = prod(p^amostra * (1-p)^(1-amostra)) return(saida) } log.verossimilhanca <- function(p,amostra){ saida <- sum(amostra*log(p)) + sum((1-amostra)*log(1-p)) return(saida) } ## Já sabemos que p esta entre 0 e 1 então p.estimado <- seq(0,1,by=0.001) ## Podemos avaliar a verossimilhança para o conjunto de p's e então escolher o que é mais provável ter gerado a amostra. vero.10 <- apply(as.matrix(p.estimado),1,verossimilhanca,amostra=amostra10) vero.30 <- apply(as.matrix(p.estimado),1,verossimilhanca,amostra=amostra30) vero.50 <- apply(as.matrix(p.estimado),1,verossimilhanca,amostra=amostra50) vero.100 <- apply(as.matrix(p.estimado),1,verossimilhanca,amostra=amostra100) vero.1000 <- apply(as.matrix(p.estimado),1,verossimilhanca,amostra=amostra1000) ## Podemos fazer um gráfico plot(p.estimado ,vero.10/max(vero.10),type="l",ylab="Verossimilhança",xlab="p",col=1,lty=1) lines(p.estimado ,vero.30/max(vero.30),type="l",ylab="Verossimilhança",xlab="p",col=2,lty=2) lines(p.estimado ,vero.50/max(vero.50),type="l",ylab="Verossimilhança",xlab="p",col=3,lty=3) lines(p.estimado ,vero.100/max(vero.100),type="l",ylab="Verossimilhança",xlab="p",col=4,lty=4) lines(p.estimado ,vero.1000/max(vero.1000),type="l",ylab="Verossimilhança",xlab="p",col=5,lty=5) legend("topright", legend = c("n=10 ", "n=30", "n=50","n=100","n=1000"), lty = 1:5, col=1:5) ## Vamos dar um zoom plot(p.estimado ,vero.10/max(vero.10),type="l",ylab="Verossimilhança",xlab="p",col=1,lty=1,xlim=c(0.2,0.6)) lines(p.estimado ,vero.30/max(vero.30),type="l",ylab="Verossimilhança",xlab="p",col=2,lty=2,xlim=c(0.2,0.6)) lines(p.estimado ,vero.50/max(vero.50),type="l",ylab="Verossimilhança",xlab="p",col=3,lty=3,xlim=c(0.2,0.6)) lines(p.estimado ,vero.100/max(vero.100),type="l",ylab="Verossimilhança",xlab="p",col=4,lty=4,xlim=c(0.2,0.6)) lines(p.estimado ,vero.1000/max(vero.1000),type="l",ylab="Verossimilhança",xlab="p",col=5,lty=5,xlim=c(0.2,0.6)) legend("topright", legend = c("n=10 ", "n=30", "n=50","n=100","n=1000"), lty = 1:5, col=1:5) ## Vamos obter o Máximo estima <- data.frame(p.estimado,n10 = vero.10/max(vero.10), n30 = vero.30/max(vero.30),n50 = vero.50/max(vero.50), n100 = vero.100/max(vero.100),n1000=vero.1000/max(vero.1000)) maxi10 <- estima[which(estima$n10 == max(estima$n10)),] maxi30 <- estima[which(estima$n30 == max(estima$n30)),] maxi50 <- estima[which(estima$n50 == max(estima$n50)),] maxi100 <- estima[which(estima$n100 == max(estima$n100)),] maxi1000 <- estima[which(estima$n1000 == max(estima$n1000)),] resultados <- rbind(maxi10,maxi30,maxi50,maxi100,maxi1000) ### Podemos obter estes valores por um algoritmo numerico (otimizador) est10 = optimize(f = verossimilhanca , interval =c(0,1) , amostra=amostra10,maximum=TRUE) est30 = optimize(f = verossimilhanca , interval =c(0,1) , amostra=amostra30,maximum=TRUE) est50 = optimize(f = verossimilhanca , interval =c(0,1) , amostra=amostra50,maximum=TRUE) est100 = optimize(f = verossimilhanca , interval =c(0,1) , amostra=amostra100,maximum=TRUE) est1000 = optimize(f = verossimilhanca , interval =c(0,1) , amostra=amostra1000,maximum=TRUE) ### Resultados obtidos pelo algoritmo numerico resultados$optim = c(est10[1],est30[1],est50[1],est100[1],est1000[1]) ## Como seria isso analiticamente resultados$analitico = c(mean(amostra10),mean(amostra30),mean(amostra50),mean(amostra100),mean(amostra1000)) resultados ## E se ao invés de uma amostra tivessemos várias amostra <- matrix(ncol=1000,nrow=100) for(i in 1:1000){ amostra[,i] <- rbinom(100,size = 1, prob=p)} medias = (apply(amostra,2,mean)) vero <- matrix(ncol=1000,nrow=1001) for(i in 1:1000){ vero[,i] <- apply(as.matrix(p.estimado),1,verossimilhanca,amostra=amostra[,i])} plot(p.estimado,vero[,1]/max(vero[,1]),type="l") for(i in 2:1000){ lines(p.estimado,vero[,i]/max(vero[,i]),col=i,lty=i) } points(x=p,y=1,pch=18) rug(jitter(medias,amount=0.008)) %abline(v=medias) ### Distribuição do estimador de Máxima Verossimilhança hist(apply(amostra,2,mean),ylim=c(0,10),prob=TRUE) curve(dnorm(x,mean=mean(medias),sd=sqrt(mean(medias)*(1-mean(medias))/100)),add=TRUE) curve(dnorm(x,mean=p,sd=sqrt(p*(1-p)/100)),add=TRUE,col="black",lty=2) lines(density(apply(amostra,2,mean)),col="blue",lty=3) quantile(medias,c(0.025,0.975)) erro.padrao = qnorm(0.975)*sqrt(mean(medias)*(1-mean(medias))/100) c(mean(medias)-erro.padrao,mean(medias),mean(medias)+erro.padrao) abline(v=0.35,col="red",lwd=3) abline(v=0.56,col="red",lwd=3) mean(apply(amostra,2,mean)) p ## Refazendo com um zoom da área mais densa do gráfico plot(p.estimado,vero[,1]/max(vero[,1]),type="l",xlim=c(0.2,0.6)) for(i in 2:1000){ lines(p.estimado,vero[,i]/max(vero[,i]),col=i,lty=i) } points(x=p,y=1,pch=18) points(y=1,x=mean(apply(amostra,2,mean)),pch=20,col="red") ############################################################################################################################### ############################### FIM DO EXEMPLO ################################################################################ ############################################################################################################################### ######################### População