O R é um software livre de larga utilização em Estatística e seu número de usuários aumenta a cada dia. O Objetivo deste texto é, de forma simples, descrever de que forma alternativa podemos realizar analises presentes no livro Estatística Básica(Bussab e Morettin, 2017) utilizando o software livre R. Como o R é de código aberto ele possui um amplo número de implementações que nos auxiliam em tal tarefa. Por conta disso, as soluções nem sempre são únicas, no sentido de que muitas vezes possuímos a mesma função implementada (não necessariamente da mesma forma), em diversos pacotes diferentes. Sempre que possível, descreveremos diferentes formas de realizar as análises, para que o leitor tenha a opção de escolher aquela que melhor lhe agrada. É importante ainda estar familiarizado com a linguagem R. Você pode encontrar detalhes sobre ela nos manuais do R. Uma série de outras referências e comandos básicos também pode ser encontrada aqui. Se o leitor desejar se aprofundar na utilização do R e dos exemplos do livro com este programa, há uma página detalhada deste conteúdo com comandos e saídas que pode ser acessada clicando aqui.
library(akima)
library(boot)
library(bootstrap)
library(colorRamps)
library(MASS)
library(ggExtra)
library(ggplot2)
library(gmodels)
library(knitr)
library(mvtnorm)
library(pander)
library(reshape2)
library(scales)
library(scatterplot3d)
library(stats)
library(stats4)
library(tables)
library(xtable)Se o leitor deseja apenas acessar os scripts das funções programadas num arquivo texto(extensão .R), clique aqui. Os dados e funções programadas em uma imagem do R (arquivo .RData) podem ser acessados aqui ou simplesmente colando a seguinte linha no console do R:
load(url(description = "https://www.ime.usp.br/~pam/dados.RData"))tab2_1<-read.table("tabela2_1.csv", dec=",", sep=";",h=T)
names(tab2_1)
summary(tab2_1$salario)ni<-table(tab2_1$grau_instrucao) # Calcula a tabela de frequências absolutas e armazena o resultado em 'mytab'
fi<-prop.table(ni) # Tabela de frequências relativas (f_i)
p_fi<-100*prop.table(ni) # Porcentagem (100 f_i)
# Adiciona linhas de total
ni<-c(ni,sum(ni))
fi<-c(fi,sum(fi))
p_fi<-c(p_fi,sum(p_fi))
names(ni)[4]<-"Total"
tab2_2<-cbind(ni,fi=round(fi,digits=2),p_fi=round(p_fi,digits=2))
tab2_2#quebras de linha apenas ilustrativas para facilitar a leitura
tab2_3<-as.data.frame(
t(rbind(
ni=c(650,1020,330,2000),
p_fi=c(32.5,51,16.5,1)
))
,row.names =c("Fundamental","Médio","Superior","Total")
)
tab2_3ni<-table(cut(tab2_1$salario, breaks = seq(4,24,by=4),right=FALSE)) # Frequencias por categorias
tab2_4 <- rbind(ni, p_fi = 100*prop.table(ni)) # Frequencias relativas em %
#quebras de linha apenas ilustrativas para facilitar a leitura
tab2_4 <- as.data.frame(
t(cbind(
tab2_4,
c(sum(tab2_4[1,]),sum(tab2_4[2,])
))),row.names =c(colnames(tab2_4),"Total")) #Construcao da tabela
tab2_4<-transform(tab2_4,p_fi=round(p_fi,digits=2))
tab2_4#quebras de linha apenas ilustrativas para facilitar a leitura
barplot(
table(tab2_1$grau_instrucao),
ylab="Frequência",
cex.names=0.7,
names.arg = c("Fundamental","Médio", "Superior"),
col="darkgrey",
border=NA,
#main="Figura 2.2: Gráfico em barras para a variável Y: grau de instrução.",
axes=TRUE,
ylim=c(0,20)
)labs<-paste(1:3,"(",tab2_2[1:3,1],";",round(tab2_2[1:3,3],1),"%)",sep="")
pie(table(tab2_1$grau_instrucao),labels=labs)
#title("Figura 2.3: Gráfico em setores para a variável Y: grau de instrução")
legend(-1.1,-0.8,legend=c("1=Fundamental, 2=Médio, 3=Superior"),border=NA,box.col=NA)#quebras de linha apenas ilustrativas para facilitar a leitura
barplot(
table(tab2_1$n_filhos),
ylab="Frequência",
cex.names=0.7,
col="darkgrey",
#main="Figura 2.4: Gráfico em barras para a variável Z: Numero de filhos.",
border=NA)par(mfrow=c(1,3),pin=c(2,2))
tbs<-as.data.frame(cbind(x=0:5,y=c(1,1,1,1,NA,1)),row.names = as.integer(c(4,5,7,3,NA,1)))
plot(tbs,ylim=c(0,7),pch=19,ylab=NA,bty="n",yaxt="n", col="darkblue",xlab="(a)")
text(x=tbs$x,y=tbs$y,rownames(tbs),pos=3)
stripchart(tab2_1$n_filhos,method = "stack", offset = 1, pch = 19, col="darkblue",ylim=c(0,7),ylab=NA,bty="n",yaxt="n",xlab="(b)",cex=1)
plot(table(tab2_1$n_filhos),type="p", col="darkblue",pch = 19,bty="n",ylab=NA,xlab="(c)")ni<-table(tab2_1$n_filhos) # Frequencias absolutas
tab2_5 <- rbind(ni, p_fi = 100*prop.table(ni)) # Frequencias relativas em %
#quebras de linha apenas ilustrativas para facilitar a leitura
tab2_5 <- as.data.frame(
t(cbind(
tab2_5,
c(sum(tab2_5[1,]),sum(tab2_5[2,])
))),row.names =c(colnames(tab2_5),"Total")) #Construcao da tabela
tab2_5<-transform(tab2_5,p_fi=round(p_fi,digits=2))
tab2_5#quebras de linha apenas ilustrativas para facilitar a leitura
barplot(
table(cut(tab2_1$salario, breaks = seq(4,24,by=4),right=FALSE)),
ylab="Frequência",
xaxt="n",
cex.names=0.7,
col="darkgrey",
border=NA)
#main="Figura 2.6: Gráfico em barras para a variável S: salários" )
axis(1,at=c(.75,1.9,3.1,4.3,5.5),labels=seq(6,22,4),tick=F)ni<-table(cut(tab2_1$salario, breaks = seq(4,24,by=4),right=FALSE)) # Frequencias por categorias
#quebras de linha apenas ilustrativas para facilitar a leitura
tab2_6 <- rbind(si=seq(6,22,by=4),ni, p_fi = 100*prop.table(ni)) # Frequencias relativas em %
tab2_6 <- as.data.frame(
t(cbind(tab2_6,
c(NA,sum(tab2_6[2,]),sum(tab2_6[3,])
))),row.names =c(colnames(tab2_6),"Total")) #Construcao da tabela
tab2_6<-transform(tab2_6,p_fi=round(p_fi,digits=2))
tab2_6fig27<-hist(tab2_1$salario, breaks = seq(4,24,by=4),right=FALSE,probability = T,plot=F)
aux<-with(fig27, 100 * density* diff(breaks)[1])
labs <- paste(round(aux), "%", sep="")
#quebras de linha apenas ilustrativas para facilitar a leitura
plot(fig27,
freq = FALSE, labels = labs,
ylab="Densidade de Frequência",
xlab="Salário",
col="darkgrey",
border="white",
#labels=T,
#main="Figura 2.7: Histograma da variável S: salários",
xlim=c(0,24), xaxp=c(0,24,6),
ylim=c(0,.1))fig28<-hist(tab2_1$n_filhos, right=F, breaks=seq(-.5,5.5,1),plot=F)
aux<-with(fig28, 100 * density* diff(breaks)[1])
labs <- paste(round(aux), "%", sep="")
#quebras de linha apenas ilustrativas para facilitar a leitura
plot(fig28,
ylab="Densidade de Frequência",
xlab="Número de Filhos",
col="darkgrey",
border="white",
bty="n",yaxt="n",ylim=c(0,8),
#main="Figura 2.8: Histograma da variável Z: número de filhos"
labels=labs)#print("Figura 2.9: Ramo-e-folhas para a Variável S: salários.")
stem(tab2_1$salario,scale=2)#quebras de linha apenas ilustrativas para facilitar a leitura
dureza<-c(53 ,70.2,84.3,69.5,77.8,87.5,53.4,82.5,67.3,54.1,
70.5,71.4,95.4,51.1,74.4,55.7,63.5,85.8,53.5,64.3,
82.7,78.5,55.7,69.1,72.3,59.5,55.3,73 ,52.4,50.7
)
stem(as.integer(dureza),scale=.5)#print("Figura 2.11: Ramo-e-folhas para dados de dureza, com ramos divididos.")
stem(as.integer(dureza),scale=1)cd_notas<-read.table("cd-notas.csv",h=T,skip=4,sep=";",dec=",")
hist(cd_notas$nota,col="darkgrey",
#main="Figura 2.12: Histograma para o CD-Notas. R.",
xlab="Notas",ylab="Frequência",border="white")Figura 2.12: Histograma para o CD-Notas. R.
#### Figura 2.13
stripchart(cd_notas$nota,method = "stack", offset = 2, at=0,
#main="Figura 2.13: Gráfico de dispersão unidimensional para CD-Notas. R.",
pch = 19, col="darkblue",ylab=NA,cex=0.5)stem(cd_notas$nota)cd_poluicao<-read.table("cd-poluicao.csv",h=T,skip=8,sep=";",dec=",")
#quebras de linha apenas ilustrativas para facilitar a leitura
plot.ts(cd_poluicao$temp,
#main="Figura 2.15: Dados de Temperatura de São Paulo. R.",
xlab="Dia", ylab="Grau",col="darkgrey")#quebras de linha apenas ilustrativas para facilitar a leitura
hist(cd_poluicao$temp,col="darkgrey",xlab="Temperatura",border="white",
#main="Figura 2.16: Histograma dos dados de Temperatura de São Paulo. R.",
ylab="")stripchart(cd_poluicao$temp, method = "stack", offset = 2, at=0,
#main="Figura 2.17: Gráfico de dispersão unidimensional \n para os dados de temperatura de São Paulo. R.",
pch = 19, col="darkblue",ylab=NA,cex=0.5)#print("Figura 2.18: Ramo-e-Folha para os dados de temperatura de São Paulo. R.")
stem(cd_poluicao$temp, scale=.5)#mean(x, na.rm = FALSE, ...)
# na.rm: Se este argumento for TRUE, então os valores NA são desconsiderados para o calculo da média.
mean(tab2_1$n_filhos,na.rm=T)fig27<-hist(tab2_1$salario, breaks = seq(4,24,by=4),right=FALSE, plot=F)
median(tab2_1$salario)
moda2(tab2_1$grau_instrucao)
tab2_1$grau_instr_num<-1*c(tab2_1$grau_instrucao=="ensino fundamental")+2*c(tab2_1$grau_instrucao=="ensino médio")+3*c(tab2_1$grau_instrucao=="superior")
median(tab2_1$grau_instr_num)
levels(as.factor(tab2_1$grau_instrucao))[median(tab2_1$grau_instr_num)]
mean(tab2_1$salario[
tab2_1$salario>quantile(tab2_1$salario,probs=0.01)
&
tab2_1$salario<quantile(tab2_1$salario,probs=0.99)
])
#Note que este valor é menor que o valor da média da variável sem a "truncagem":
mean(tab2_1$salario)n<-sum(!is.na(tab2_1$n_filhos))
#Variância:
var(tab2_1$n_filhos,na.rm=T)*(n-1)/n
#Desvio padrão:
sd(tab2_1$n_filhos,na.rm=T)*(n-1)/n
#Desvio médio:
mean(abs(tab2_1$n_filhos[!is.na(tab2_1$n_filhos)]-mean(tab2_1$n_filhos,na.rm=T)))
# Exemplo:
x<-c(0,1,NA,NA,0,-2,pi,NA)
is.na(x)
!is.na(x)mean(tab2_1$salario,na.rm=T)
dm(tab2_1$salario)
#n<-sum(!is.na(tab2_1$n_filhos))
varp(tab2_1$salario,na.rm=T)x<-c(15,5,3,8,10,2,7,11,12)
median(x)
mean(x)
quantile(x,probs=c(.25,.50,.75),type=5) # "type 6" indica o mesmo método de cálculo que o SPSS e Minitab
#Agora, adicionando a observação 67 ao vetor x, temos:
x<-c(sort(x),67) # A função `sort` ordena o vetor
x
median(x)
mean(x)
quantile(x,probs=c(.25,.50,.75),type=5)
# 20% das observações a esquerda
quantile(x,probs=0.2)fig27<-hist(tab2_1$salario, breaks = seq(4,24,by=4),right=FALSE,probability = T,plot=F)
aux<-with(fig27, 100 * density* diff(breaks)[1])
labs <- paste(round(aux), "%", sep="")
#quebras de linha apenas ilustrativas para facilitar a leitura
plot(fig27,
freq = FALSE, labels = labs,
xlab="Salário",
ylab="",
col="darkgrey",
border="white",
yaxt="n",
xlim=c(0,24), xaxp=c(0,24,6),
ylim=c(0,.1), main="")
median(tab2_1$salario)
quantile(tab2_1$salario,probs=c(.20,.5,.95,.75))
q<-quantile(tab2_1$salario,probs=c(.25,.5,.75)) # armazenando as informações de quantis em um objeto `q`
q
dq=q[3]-q[1]
print(paste("Distância inter-quartis (d) = ", dq))# Utilizando os dados do exemplo 3.5 :
x<-c(15,5,3,8,10,2,7,11,12)
summary(x)cd_municipios<-read.table("cd-municipios.csv",h=T,skip=4,sep=";",dec=",")
# Medidas resumo para a base de dados dos municípios (n=30)
summary(cd_municipios$populacao,digits = 4)
# Medidas resumo para a base de dados dos municípios, desconsiderando São Paulo e Rio de Janeiro (n=28)
summary(cd_municipios$populacao[-c(1,2)],digits = 4) munic<-cd_municipios[order(cd_municipios$populacao,decreasing = TRUE),]
munic
munic<-cd_municipios[1:15,]
munic
# Medidas resumo para a base de dados dos 15 municípios mais populosos (n=15)
summary(munic$populacao)boxplot(munic$populacao,
pch="*", # tipo de marcador dos outliers
col="lightblue", # cor do preenchimento do box plot
border="darkgrey", # cor da linha do box plot
boxwex=0.3 # Tamanho da caixa
)
# Os comandos 'text'a seguir imprimem no box plot os nomes das cidades dos 4 pontos destacados
text(x=cd_municipios$populacao[1],label=cd_municipios$municipio[1],pos=4,cex=0.7) # Máximo
text(x=cd_municipios$populacao[2],label=cd_municipios$municipio[2],pos=4,cex=0.7) # 2a. maior observacao
text(x=cd_municipios$populacao[3],label=cd_municipios$municipio[3],pos=4,cex=0.7) # 3a. maior observacao
text(x=cd_municipios$populacao[15],label=cd_municipios$municipio[15],pos=4,cex=0.7) # 15a. observacao
#Para imprimir o box plot desconsiderando aqueles pontos que foram considerados outliers de acordo com seu critério, faça:
boxplot(munic$populacao, col="lightblue", border="darkgrey", boxwex=0.3, outline=FALSE)
title("Box plot para os quinze maiores \n municípios do Brasil, desconsiderando outliers")x<-c(0.5,2.3,4,6.4,8,15.3,13.5,12,9.8)
# Mediana de x
median(x)
# Calculando os ui e vi
u<-sort((median(x)-x)[(median(x)-x)>0],decreasing = T)
v<-sort((x-median(x))[(x-median(x))>0],decreasing = T)
u
vu<-median(cd_municipios$populacao)-cd_municipios$populacao
v<-cd_municipios$populacao-median(cd_municipios$populacao)
plot(sort(u),sort(v), pch=19, xlab="ui", ylab="vi",col="darkblue",xlim=c(0,max(u)),ylim=c(0,max(v)))
#title("Figura 3.11: Gráfico de simetria para o CD-Municípios.")
abline(0,1)xp<-list() # declarando xp como uma lista vazia
xp[[1]]<-log(cd_municipios$populacao) # Transformação log
xp[[2]]<-(cd_municipios$populacao)^(1/4) # Transformação p=1/4
xp[[3]]<-(cd_municipios$populacao)^(1/2) # Transformação p=1/2
xp[[4]]<-(cd_municipios$populacao)^(1/3) # Transformação p=1/3
par(mfrow=c(2,2))
hist(xp[[1]], main="log", ylab="", xlab="", col="darkgrey", border="white",cex.axis=0.8,cex.main=0.8)
hist(xp[[2]], main="p=1/4", ylab="", xlab="", col="darkgrey", border="white",cex.axis=0.8,cex.main=0.8)
hist(xp[[3]], main="p=1/2", ylab="", xlab="", col="darkgrey", border="white",cex.axis=0.8,cex.main=0.8)
hist(xp[[4]], main="p=1/3", ylab="", xlab="", col="darkgrey", border="white",cex.axis=0.8,cex.main=0.8)boxplot(xp,pch="-", col="lightblue", border="darkgrey", boxwex=0.3, names=c("p=0","p=1/4","p=1/2","p=1/3"))#Podemos construir uma nova função, digamos summary2, que conterá mais informação do que a summary atual:
summary2<-function(x,limtrim=c(.01,.99), pop=FALSE, na.rm=TRUE, ...){
# x : argumento numerico para o qual calcularemos as estatísticas descritivas
# limtrim: Limites para a média truncada
# pop: se esta variável é TRUE, entao utiliza-se 'n' como denominador da estimativa da variância. caso contrário, 'n-1'.
# na.rm: ação para os outliers
pop=TRUE
na.rm=TRUE
newsumm<-matrix(NA,nrow=dim(as.matrix(x))[2],ncol=11)
x<-as.matrix(x)
newsumm[,1]<-apply(x,2,function(...){return(dim(x)[1])})
newsumm[,2:7]<-t(apply(as.data.frame(x),2,summary))
newsumm[,8]<-apply(x,2,function(...){mean(x[x>quantile(x,probs=0.01) &
x<quantile(x,probs=0.99)])})
if (pop==TRUE & na.rm==TRUE){ # calcula var_n desconsiderando os missings
n<-sum(!is.na(x))
newsumm[,9]<-var(x,na.rm=na.rm)*(n-1)/n
}
else
{
if(pop==TRUE & na.rm==FALSE ){
n<-length(x)
newsumm[,9]<-var(x)*(n-1)/n
}
else {
if(pop==FALSE & na.rm==TRUE){
n<-sum(!is.na(x))
newsumm[,9]<-var(x,na.rm=TRUE)
}
else
{
n<-length(x)
newsumm[,9]<-var(x)
}
}
}
newsumm[,1]<-n
newsumm[,10]<-sqrt(newsumm[,9])
newsumm[,11]<-sqrt(newsumm[,10]/newsumm[,1])
colnames(newsumm)<-c("N", "Min.","1st Qu.","Median","Mean","3rd Qu.","Max.","Tr Mean","Var","StDev","SE Mean")
return(t(newsumm))
}
cd_notas<-read.table("cd-notas.csv",header=TRUE,skip=4,sep=";",dec=",") # Leitura dos dados
summary2(cd_notas$nota)boxplot(cd_notas$nota,pch="-", col="lightblue", border="darkgrey")u<-median(cd_notas$nota)-cd_notas$nota
v<-cd_notas$nota-median(cd_notas$nota)
plot(sort(u),sort(v), pch=19, xlab="ui", ylab="vi",col="darkblue",xlim=c(0,max(u)),ylim=c(0,max(v)))
title("Figura 3.15: Gráfico de simetria para o CD-Notas.")
abline(0,1)# Quadro 3.4 Medidas descritivas para temperaturas. R.
summary(cd_poluicao$temp)boxplot(cd_poluicao$temp,pch="-", col="lightblue", border="darkgrey")u<-median(cd_poluicao$temp)-cd_poluicao$temp
v<-cd_poluicao$temp-median(cd_poluicao$temp)
plot(sort(u),sort(v), pch=19, xlab="ui", ylab="vi",col="darkblue",xlim=c(0,max(u)),ylim=c(0,max(v)))
abline(0,1)x<-c(15,5,3,8,10,2,7,11,12)
plot(ecdf(x),xlim=c(0,16),ylim=c(0,1),main="",ylab="",pch=20,cex=0.7) # 'ecdf' é a função que calcula a função distribuição empírica
lines(y = c(0.5,0.5),x=c(-1,8),col="grey")
lines(y = c(0,0.5),x=c(8,8),col="grey")
text(-.5,0.5,"p=9/18",pos=4,col="darkblue")
text(8,0,"8",pos=3,col="darkblue")
#title("Figura 3.18: Funções de distribuição empírica(F_e) \n e f.d.e. alisada (Fs_e) para o Exemplo 3.5")
par(new=T) # Utilize este comando para sobrepor o gráfico a seguir no gráfico anterior
# Construiremos a seguir uma função 'scdf', que será a ecdf acima suavisada
scdf<-function(x){
f<-t<-table(x)
st<-sum(t[1:length(x)])
for(i in 1:length(x))
f[i]<-sum(t[1:i])/st - (t[i]/st)/2
return(f)
}
par(new=T)
plot(scdf(x),pch=20,cex=0.7,lty=2, type="b",xlab="",ylab="",col="grey",xaxt="n",yaxt="n",xlim=c(0,16),ylim=c(0,1))
legend(12,.4,c("F_e","Fs_e"),lty=c(1,2),col=c("black","grey"),lwd=2,cex=0.8)attach(tab2_1)
tab4_2<-table(reg_procedencia,grau_instrucao)
tab4_2
Total_linha<-margin.table(tab4_2,2) # O argumento 2 define a marginal da linha
Total_coluna<-margin.table(tab4_2,1) # O argumento 1 define a marginal da coluna
tab4_2_<-rbind(cbind(tab4_2,Total_coluna),c(Total_linha, sum(Total_coluna)))
dimnames(tab4_2_)[[1]][4]<-"Total_linha"
tab4_2_
Total_linha<-apply(tab4_2, MARGIN=2, FUN=sum)
Total_coluna<-apply(tab4_2, MARGIN=1, FUN=sum)
tab4_2_<-rbind(cbind(tab4_2,Total_coluna),c(Total_linha, sum(Total_coluna)))
dimnames(tab4_2_)[[1]][4]<-"Total_linha" tab4_3<-prop.table(tab4_2)
tab4_3
Total_linha<-margin.table(tab4_3,2)
Total_coluna<-margin.table(tab4_3,1)
tab4_3_<-rbind(cbind(tab4_3,Total_coluna),c(Total_linha, sum(Total_coluna)))
dimnames(tab4_3_)[[1]][4]<-"Total_linha"
tab4_3_tab4_4<-prop.table(tab4_2,2)
tab4_4
Total_linha<-margin.table(tab4_4,2)
Total_coluna<-margin.table(tab4_2,1)/sum(margin.table(tab4_2,1))
tab4_4_<-rbind(cbind(tab4_4,Total_coluna),c(Total_linha, sum(Total_coluna)))
dimnames(tab4_4_)[[1]][4]<-"Total_linha"
tab4_4_gmodels#Tabela 4.2
CrossTable(reg_procedencia,grau_instrucao,
prop.r=FALSE, # Se TRUE, entao retorna as proporções nas linhas
prop.c=FALSE, # Se TRUE, entao retorna as proporções nas colunas
prop.t=FALSE, # Se TRUE, entao retorna as proporções em relação ao total
prop.chisq=FALSE # Se TRUE, entao retorna a contribuição de cada casela para a estatística de Qui-quadrado
)
#Tabela 4.3
CrossTable(reg_procedencia,grau_instrucao,prop.r=FALSE, prop.c=FALSE,
prop.t=TRUE, prop.chisq=FALSE)
#Tabela 4.4
CrossTable(reg_procedencia,grau_instrucao,prop.r=FALSE, prop.c=TRUE,
prop.t=FALSE, prop.chisq=FALSE)ggplot(melt(tab4_2_[1:3,],value.name = "contagem",
varnames = c("reg_procedencia","grau_instrucao") ))+ ## `melt` empilha os dados no formato necessário para o ggplot
aes(x=grau_instrucao,y=contagem,fill=reg_procedencia) + ## Variáveis a serem plotadas.
geom_bar(stat="identity", position = "fill") + ## Define o gráfico de barras percentual empilhado
scale_fill_brewer(name="Região de\n Procedência")+ ## Opções do preenchimento do gráfico (label e paleta de cores)
scale_y_continuous(labels = percent_format()) + ## Formato do eixo Y em porcentagem
theme_bw()+ ## Define a cor do fundo do gráfico: neste caso, branco
#theme(legend.position="bottom") + ## Define a posição da legenda abaixo do gráfico
#ggtitle("Figura 4.1: Distribuição da região de procedência por grau de instrução")+
xlab("Grau de Instrução") + ylab("") ## Define os `labels` dos eixostab4_5<-as.table(matrix(c(85,55,35,25),ncol=2))
dimnames(tab4_5)[[1]] = c("Economia", "Administração")
dimnames(tab4_5)[[2]] = c("Masculino", "Feminino")
tab4_5
Total_linha<-margin.table(tab4_5,2)
Total_coluna<-margin.table(tab4_5,1)
tab4_5_<-rbind(cbind(tab4_5,Total_coluna),c(Total_linha, sum(Total_coluna)))
dimnames(tab4_5_)[[1]][3]<-"Total_linha" tab4_6<-prop.table(tab4_5,2)
tab4_6
Total_linha<-margin.table(tab4_6,2)
Total_coluna<-margin.table(tab4_6,1)/sum(margin.table(tab4_6,1))
tab4_6_<-rbind(cbind(tab4_6,Total_coluna),c(Total_linha, sum(Total_coluna)))
dimnames(tab4_6_)[[1]][3]<-"Total_linha" # Construindo os dados do exemplo:
dados.tab4_7<-data.frame(rbind(
matrix(rep(c("1.Masculino","1.Física"),times=100),ncol=2,byrow=T),
matrix(rep(c("2.Feminino","1.Física"),times=20),ncol=2,byrow=T),
matrix(rep(c("1.Masculino","2.Ciências Sociais"),times=40),ncol=2,byrow=T),
matrix(rep(c("2.Feminino","2.Ciências Sociais"),times=40),ncol=2,byrow=T)))
colnames(dados.tab4_7)<-c("sexo","curso")
CrossTable(dados.tab4_7$curso,dados.tab4_7$sexo,prop.r=FALSE, prop.c=TRUE,
prop.t=FALSE, prop.chisq=FALSE)# Construindo os dados do exemplo:
dados.tab4_8<-data.frame(rbind(
matrix(rep(c("1.Consumidor","1.São Paulo"),times=214),ncol=2,byrow=T),
matrix(rep(c("1.Consumidor","2.Paraná"),times=51),ncol=2,byrow=T),
matrix(rep(c("1.Consumidor","3.Rio G. do Sul"),times=111),ncol=2,byrow=T),
matrix(rep(c("2.Produtor","1.São Paulo"),times=237),ncol=2,byrow=T),
matrix(rep(c("2.Produtor","2.Paraná"),times=102),ncol=2,byrow=T),
matrix(rep(c("2.Produtor","3.Rio G. do Sul"),times=304),ncol=2,byrow=T),
matrix(rep(c("3.Escola","1.São Paulo"),times=78),ncol=2,byrow=T),
matrix(rep(c("3.Escola","2.Paraná"),times=126),ncol=2,byrow=T),
matrix(rep(c("3.Escola","3.Rio G. do Sul"),times=139),ncol=2,byrow=T),
matrix(rep(c("4.Outras","1.São Paulo"),times=119),ncol=2,byrow=T),
matrix(rep(c("4.Outras","2.Paraná"),times=22),ncol=2,byrow=T),
matrix(rep(c("4.Outras","3.Rio G. do Sul"),times=48),ncol=2,byrow=T)))
colnames(dados.tab4_8)<-c("tipo_de_cooperativa","estado")
attach(dados.tab4_8)CrossTable(estado,tipo_de_cooperativa,
prop.r=TRUE, prop.c=FALSE, prop.t=FALSE, prop.chisq=FALSE,
digits=2)CrossTable(estado,tipo_de_cooperativa,
prop.r=FALSE, prop.c=FALSE, prop.t=FALSE, prop.chisq=FALSE, expected=TRUE,
digits=0)CrossTable(estado,tipo_de_cooperativa,
prop.r=FALSE, prop.c=FALSE, prop.t=FALSE, resid=TRUE, prop.chisq=TRUE,
format="SPSS", digits=2)
tab4_8<-table(estado,tipo_de_cooperativa)
testequi<-chisq.test(tab4_8)
testequi### (4.5) ###
C = sqrt(testequi$statistic/(testequi$statistic+1551))
round(C,digits=2) # Arredondando para duas casas decimais
### (4.6) ###
T = sqrt(testequi$statistic/(1551*(4-1)*(3-1)))
round(T,digits=2) # Arredondando para duas casas decimaisdados.ex4_4<-data.frame(agente=c("A","B","C","D","E","F","G","H","I","J"),
anos_servico=c(2,3,4,5,4,6,7,8,8,10),
n_clientes=c(48,50,56,52,43,60,62,58,64,72))print(kable(dados.ex4_4, caption="**Tabela 4.12** Número de anos de serviço(X) por número de clientes(Y) de agentes de uma companhia de seguros."))attach(dados.ex4_4)
plot(n_clientes,anos_servico,pch=20, col="darkblue",
xlab="Anos de Serviço", ylab="Número de Clientes",cex=2)cd_brasil<-read.table("cd-brasil.csv",h=T,skip=7,sep=";")
names(cd_brasil)
attach(cd_brasil)
plot(pop_urbana,pop_rural,pch=20, col="darkblue",
xlab="População urbana",ylab="População rural",cex=2)dados.ex4_6<-data.frame(
familia=c("A","B","C","D","E","F","G","H","I","J"),
renda_bruta=c(12,16,18,20,28,30,40,48,50,54),
gasto_saude=c(7.2,7.4,7,6.5,6.6,6.7,6,5.6,6,5.5))
attach(dados.ex4_6)
plot(renda_bruta,gasto_saude,pch=20, col="darkblue",
xlab="Renda Bruta",ylab="% de gasto com saúde",cex=2)dados.tab4_14<-data.frame(
individuo=c("A","B","C","D","E","F","G","H"),
resultado=c(45,52,61,70,74,76,80,90),
tempo=c(343,368,355,334,337,381,345,375))
attach(dados.tab4_14)
plot(resultado,tempo,pch=20, col="darkblue",
xlab="Resultado teste",ylab="Tempo",cex=2)attach(dados.ex4_4)
tab4_15<-dados.ex4_4
tab4_15$dsvx<-anos_servico-mean(anos_servico) # Desvio da variável X em relação a sua média
tab4_15$dsvy<-n_clientes-mean(n_clientes) # Desvio da variável Y em relação a sua média
tab4_15$zx<-tab4_15$dsvx/sqrt(varp(anos_servico)) # Desvios de X padronizados pelo seu desvio-padrão
tab4_15$zy<-tab4_15$dsvy/sqrt(varp(n_clientes)) # Desvios de Y padronizados pelo seu desvio-padrão
tab4_15$zxzy<-tab4_15$zx*tab4_15$zy # Coeficiente de correlação entre X e Y
attach(tab4_15)
par(mfrow=c(1,2))
plot(dsvx,dsvy,pch=20, col="darkblue",
xlab=expression(x-bar(x)),ylab=expression(y-bar(y)))
abline(v=0, lty=2,col="grey")
abline(h=0, lty=2,col="grey")
plot(zx,zy,pch=20, col="darkblue",
xlab=expression(z[x]),ylab=expression(z[y]))
abline(v=0, lty=2,col="grey")
abline(h=0, lty=2,col="grey")
#cor(X,Y) # Calcula a correlação linear entre X e Y
#cov(X,Y) # Calcula a covariância entre X e Yattach(tab2_1)
tapply(salario,grau_instrucao,summary)-> summary.tab4_16 # Calculando medidas de posição por grupo
tapply(salario,grau_instrucao,varp)-> varp.tab4_16 # Calculando variancia por grupo
tapply(salario,grau_instrucao,length)-> n.tab4_16 # Calculando n por grupo
summary.tab4_16[[4]]<-summary(salario) # Calculando medidas de posição de toda a amostra
varp.tab4_16[[4]]<-varp(salario) # Calculando variancia de toda a amostra
n.tab4_16[[4]]<-length(salario) # Calculando n
tab4_16 <- matrix(unlist(summary.tab4_16), ncol = 6, byrow = TRUE)
tab4_16<-cbind(n.tab4_16,tab4_16,sqrt(varp.tab4_16),varp.tab4_16)
dimnames(tab4_16)<-list(
c("Fundamental","Medio","Superior", "Todos"),
c("n", "Min.","Q1","Mediana","Média","Q3","Max.","dp(S)","Var(S)"))print(kable(tab4_16,digits = 2,caption="**Tabela 4.16:** Medidas-resumo para a variável salário, segundo grau de instrução, na Companhia MB."))ggplot(tab2_1, aes(grau_instrucao, salario)) + geom_boxplot(fill = "darkblue", colour = "grey")
attach(tab2_1)
tapply(salario,reg_procedencia,summary)-> summary.tab4_17 # Calculando medidas de posição por grupo
tapply(salario,reg_procedencia,varp)-> varp.tab4_17 # Calculando variancia por grupo
tapply(salario,reg_procedencia,length)-> n.tab4_17 # Calculando n por grupo
summary.tab4_17[[4]]<-summary(salario) # Calculando medidas de posição de toda a amostra
varp.tab4_17[[4]]<-varp(salario) # Calculando variancia de toda a amostra
n.tab4_17[[4]]<-length(salario) # Calculando n
tab4_17 <- matrix(unlist(summary.tab4_17), ncol = 6, byrow = TRUE)
tab4_17<-cbind(n.tab4_17,tab4_17,sqrt(varp.tab4_17),varp.tab4_17)
dimnames(tab4_17)<-list(
c("Capital","Interior","Outra", "Todos"),
c("n", "Min.","Q1","Mediana","Média","Q3","Max.","dp(S)","Var(S)"))print(kable(tab4_17,digits = 2,caption="**Tabela 4.17:** Medidas-resumo para a variável salário segundo a região de procedência, na Companhia MB."))ggplot(tab2_1, aes(reg_procedencia, salario)) + geom_boxplot(fill = "darkblue", colour = "grey") # + geom_jitter()#R^2 para Salário por grau de instrução
varS1=t(tab4_16[1:3,1])%*%tab4_16[1:3,9]/sum(tab4_16[1:3,1])
R2.1=1-varS1/tab4_16[4,9]
print(paste("media(var(S)) :",round(varS1,2)))
print(paste("R^2 :",round(R2.1,2)))
#R^2 para Salário por região de procedência
varS2=t(tab4_17[1:3,1])%*%tab4_17[1:3,9]/sum(tab4_17[1:3,1])
R2.2=1-varS2/tab4_17[4,9]
print(paste("media(var(S)) :",round(varS2,2)))
print(paste("R^2 :",round(R2.2,2)))tab4_18<-data.frame(aluno=1:20,prova1=c(8.5,3.5,7.2,5.5,9.5,7,4.8,6.6,2.5,7,7.4,5.6,6.3,3,8.1,3.8,6.8,10,4.5,5.9),prova2=c(8,2.8,6.5,6.2,9,7.5,5.2,7.2,4,6.8,6.5,5,6.5,3,9,4,5.5,10,5.5,5))
plot(tab4_18$prova1,tab4_18$prova2,xlab="Quantis da 1a. prova",ylab="Quantis da 2a. prova",pch=16,col="darkblue")
abline(a=0,b=1)print(kable(tab4_18,caption="**Tabela 4.18:** Notas de 20 alunos em duas provas de Estatística."))cd_temperaturas<-read.table("cd-temperaturas.csv",h=T,skip=4,sep=";",dec=",") # Leitura dos dados
attach(cd_temperaturas)
q_cananeia<-quantile(cananeia,probs = seq(0,1,1/120)) # seq(0,1,1/120) constroi uma sequencia de 0 a 1 com saltos de 1/120 de distância.
q_ubatuba<-quantile(ubatuba,probs = seq(0,1,1/120))
plot(q_cananeia,q_ubatuba,xlab="Quantis Cananéia",ylab="Quantis Ubatuba",pch=16,col="darkblue")
abline(a=0,b=1)tab2_1$idade<-tab2_1$idade_anos*12+tab2_1$idade_meses
attach(tab2_1)
par(mfrow=c(1,3),pin=c(2,2))
plot(idade[grau_instrucao=="ensino fundamental"],salario[grau_instrucao=="ensino fundamental"],main="Fundamental",xlab="Idade",ylab="Salário",pch=16,col="darkblue")
plot(idade[grau_instrucao=="ensino médio"],salario[grau_instrucao=="ensino médio"],main="Médio",xlab="Idade",ylab="Salário",pch=16,col="darkblue")
plot(idade[grau_instrucao=="superior"],salario[grau_instrucao=="superior"],main="Superior",xlab="Idade",ylab="Salário",pch=16,col="darkblue")Leitura dos dados:
cd_mercado <- read.table("cd-mercado.csv",h=T,skip=4,sep=";",dec=",") # Leitura dos dados
attach(cd_mercado)
plot(telebras[1:39],indice[1:39],xlab="Telebrás",ylab="Ibovespa",pch=16,col="darkblue")
abline(lm(indice[1:39]~telebras[1:39]))
print(paste("Corr(X,Y) = ",round(cor(indice[1:39],telebras[1:39]),2)))cd_veiculos <- read.table("cd-veiculos.csv",h=T,skip=4,sep=";",dec=",") # Leitura dos dados
attach(cd_veiculos)
ggplot(cd_veiculos, aes(comprimento,preco)) + geom_point(aes(shape =N_I,colour=N_I), size = 4)A<-matrix(0,7,7) # Matriz de setas: colunas são as origens e linhas, os destinos
A[2,1]<-"2/5"
A[3,1]<-"3/5"
A[4,2]<-"2/4"
A[5,2]<-"2/4"
A[6,3]<-"3/4"
A[7,3]<-"1/4"
col <- A
col[] <- "blue"
pos=rbind(c(0.05,0.5),
c(0.3,0.33),
c(0.3,0.67),
c(0.7,0.2),
c(0.7,0.4),
c(0.7,0.6),
c(0.7,0.8)) # Matriz de coordenadas dos nós
knots=c(" ","V","B","V","B","V","B") # labels dos nós.
pp <- plotmat(A, pos = pos, name = knots,
box.cex = c(2,2,2,1,1,1,1), box.type = "none", box.prop = 1.0, # Parâmetros do nó: tamanho da letra, tipo de borda e proporção
cex.txt = 1, dtext=0, #parâmetros do texto das setas: tamanho e posição
arr.lcol = col, arr.col = col, arr.pos = 0.7, arr.type = "simple", arr.len = 0.3,
curve = 0, segment.from = 0.8, segment.to = 0.1,lwd = 2, # Parametros complementares de cor, posição da ponta da seta, tipo de seta, espessura da linha, curvatura, comprimento da linha
main = "Diagrama de árvore" # titulo
)A<-matrix(0,7,7) # Matriz de setas: colunas são as origens e linhas, os destinos
A[2,1]<-"1/2"
A[3,1]<-"1/2"
A[4,2]<-"1/2"
A[5,2]<-"1/2"
A[6,3]<-"1/2"
A[7,3]<-"1/2"
col <- A
col[] <- "blue"
pos=rbind(c(0.05,0.5),
c(0.3,0.33),
c(0.3,0.67),
c(0.7,0.2),
c(0.7,0.4),
c(0.7,0.6),
c(0.7,0.8)) # Matriz de coordenadas dos nós
knots=c(" ","R","C","R","C","R","C") # labels dos nós.
pp <- plotmat(A, pos = pos, name = knots,
box.cex = c(2,2,2,1,1,1,1), box.type = "none", box.prop = 1.0, # Parâmetros do nó: tamanho da letra, tipo de borda e proporção
cex.txt = 1, dtext=0, #parâmetros do texto das setas: tamanho e posição
arr.lcol = col, arr.col = col, arr.pos = 0.7, arr.type = "simple", arr.len = 0.3,
curve = 0, segment.from = 0.8, segment.to = 0.1,lwd = 2, # Parametros complementares de cor, posição da ponta da seta, tipo de seta, espessura da linha, curvatura, comprimento da linha
main = "Diagrama de árvore" # titulo
)lucro<-data.frame(x=c(15,10,5,-5), px=c(0.56 , 0.23, 0.02, 0.19))
plot(lucro$x,lucro$px, pch=16, col="darkblue", xlab="Lucro", ylab="p(x)", bty="n", ylim=c(0,0.7))
# Alternativamente, pode-se modificar a posição do eixo Y utilizando o comando:
# axis(2,at=seq(0,0.6,0.1),labels=seq(0,0.6,0.1),pos=0)
# Se desejar isso, é necessário adicionar a opção `yaxt="n"` no comando que gera o gráfico.o.lucro<-lucro[order(lucro$x),] # Ordenando o dataset
o.lucro<-cbind(o.lucro,fda_x=cumsum(o.lucro$px)) # Calculando a fda de X
plot(stepfun(o.lucro$x,c(0,o.lucro$fda_x)),
main="f.d.a para a v.a. X = Lucro por montagem",
xlab="Lucro", ylab="p(x)", bty="n",col="darkblue",pch=20)
lines(x=c(5,5),y=c(0,o.lucro$fda_x[1]),lty=2,col="gray")
lines(x=c(10,10),y=c(0,o.lucro$fda_x[2]),lty=2,col="gray")
lines(x=c(15,15),y=c(0,o.lucro$fda_x[3]),lty=2,col="gray")x = c(0,1,2,3)
px<-dbinom(x=x,size=3,p=1/2)
table612<-data.frame(x,px)
print(kable(table612,digits=3, align = c("r","c"),
caption="**Tabela 6.12:** Probabilidades binomiais para n=3 e p=1/2"))plot(table612$x,table612$px, pch=16, col="darkblue", xlab="Lucro", ylab="p(x)", bty="n") # x: numero de peças defeituosas desejadas
# m: Número total de peças defeituosas
# n: Número total de peças não defeituosas
# k: Número de bolas retiradas
dhyper(x=0 , m = 10, n = 90, k = 5)x<-0:9
p.5=0.5
p.005=0.005
n=10
pbinomial.5<-dbinom(x,n,p.5)
ppoisson.5<-dpois(x,lambda=n*p.5)
pbinomial.005<-dbinom(x,n,p.005)
ppoisson.005<-dpois(x,lambda=n*p.005)
table613a<-data.frame(x,pbinomial.5,ppoisson.5,pbinomial.005,ppoisson.005)
table613a<-rbind(table613a,c(10,1-sum(pbinomial.5),1-sum(ppoisson.5),1-sum(pbinomial.005),1-sum(ppoisson.005)))
print(kable(table613a))
attach(table613a)
par(mfrow=c(1,2))
plot.ts(table613a[,2:3], plot.type = "single", ylab="p(x)", xlab="x",col=c("blue","red"), main="p=0.5")
plot.ts(table613a[1:3,4:5], plot.type = "single", ylab="p(x)", xlab="x",col=c("blue","red"), main="p=0.005")x<-0:12
px<-dbinom(x,size=14,p=0.3)
fdax<-cumsum(px)
quadro61<-data.frame(x,px,fdax)
print(kable(quadro61,caption="**Quadro 6.1**: Probabilidades Binomiais geradas pelo R"))x<-0:17
px<-dpois(x,lambda=5.2)
fdax<-cumsum(px)
quadro62<-data.frame(x,px,fdax)
print(kable(quadro62,caption="**Quadro 6.1**: Probabilidades de Poisson geradas pelo R"))x<-seq(-3,3,0.1) # Cria um intervalo de -3 a 3
fdnorm<-dnorm(x = x, mean = 0, sd=1) # Calcula a fdp da distr. normal para o intervalo x
fdanorm<-pnorm(q = x, mean = 0, sd=1) # Calcula a fda da distr. normal para o intervalo x
## Imprimindo os gráficos da fdp e fda:
par(mfrow=c(1,2))
plot(x=x,y=fdnorm,type="l", col="blue",lwd=2, main="f.d.p. da Distrib. Normal padrão",xlab="z")
plot(x=x,y=fdanorm,type="l", col="blue",lwd=2, main="f.d.a. da Distrib. Normal padrão",xlab="z")
lines(x=c(0,0),y=c(0,fdanorm[x==0]),lty=2, col="gray")
par(mfrow=c(1,2))
curve(dnorm(x = x, mean = 0, sd=1), xlim=c(-3,3), col="blue",lwd=2, main="f.d.p. da Distrib. Normal padrão",xlab="z")
curve(pnorm(q = x, mean = 0, sd=1), xlim=c(-3,3), col="blue",lwd=2, main="f.d.a. da Distrib. Normal padrão",xlab="z")
lines(x=c(0,0), y=c(0,fdanorm[x==0]), lty=2, col="gray")z_a=-1.96
z_b=1.96
#Cálculo da f.d.a até z_a = -1.96
pnorm(q = z_a, mean = 0, sd = 1)
#Cálculo da f.d.a até z_b = 1.96
pnorm(q = z_b, mean = 0, sd = 1)
#Cálculo da f.d.a entre z_a = -1.96 e z_b = 1.96
pnorm(q = z_b, mean = 0, sd = 1) - pnorm(q = z_a, mean = 0, sd = 1)
par(mfrow=c(1,3))
regiao=seq(-3,-1.96,0.01)
cord.x <- c(min(regiao),regiao,max(regiao))
cord.y <- c(0,dnorm(regiao),0)
curve(dnorm(x,0,1),xlim=c(-3,3),main='Normal Padrão(até -1.96)',xlab="z")
polygon(cord.x,cord.y,col='lightgray')
regiao=seq(-3,1.96,0.01)
cord.x <- c(min(regiao),regiao,max(regiao))
cord.y <- c(0,dnorm(regiao),0)
curve(dnorm(x,0,1),xlim=c(-3,3),main='Normal Padrão(até +1.96)',xlab="z")
polygon(cord.x,cord.y,col='lightgray')
regiao=seq(-1.96,1.96,0.01)
cord.x <- c(min(regiao),regiao,max(regiao))
cord.y <- c(0,dnorm(regiao),0)
curve(dnorm(x,0,1),xlim=c(-3,3),main='Normal Padrão(de -1.96 a 1.96)',xlab="z")
polygon(cord.x,cord.y,col='lightgray')z<-seq(-3,3,0.1) # Cria um intervalo de -3 a 3
fdnorm<-dnorm(x = z, mean = 0, sd=1) # Calcula a fdp da distr. normal para o intervalo x
plot(x=z,y=fdnorm,type="l", col="blue",lwd=2, main="f.d.p. da Distrib. Normal padrão")
lines(x=c(1,1),y=c(0,fdnorm[x==1]),lty=2, col="gray") # Linha cinza tracejada em z=-1
lines(x=c(-1,-1),y=c(0,fdnorm[x==-1]),lty=2, col="gray") # Linha cinza tracejada em z=1regiao=seq(-1,3,0.01) # Regiao a ser sombreada na figura
cord.x <- c(min(regiao),regiao,max(regiao))
cord.y <- c(0,dnorm(regiao,mean=2),0)
curve(dnorm(x,2,1),xlim=c(-1,5),bty="n", main='',xlab="",ylab="",xaxt="n") #A opção `xaxt="n"` suprime o eixo x
axis(1,at=3,labels = "y") # Comando para imprimir o eixo x com labels especificos
polygon(cord.x,cord.y,col='lightgray') # imprime a regiao sombreadacurve(pnorm(q = x, mean = 0, sd=1), xlim=c(-3,3), col="blue",lwd=2,ylab="",xlab="z")
lines(x=c(0,0), y=c(0,fdanorm[x==0]), lty=2, col="gray")regiao=seq(1,2,0.01)
cord.x <- c(min(regiao),regiao,max(regiao))
cord.y <- c(0,dnorm(regiao),0)
curve(dnorm(x,0,1),xlim=c(-3,3),main='',xaxt="n",ylab="")
axis(1,at=c(1,2),labels = c("a","b"))
polygon(cord.x,cord.y,col='lightgray')regiao=seq(0,2,0.01)
cord.x <- c(min(regiao),regiao,max(regiao))
cord.y <- c(0,dnorm(regiao),0)
curve(dnorm(x,0,1),xlim=c(-3,3),main='',xaxt="n",xlab="z") #A opção `xaxt="n"` suprime o eixo x
axis(1,at=c(0,2),labels = c("0","z_c")) # Comando para imprimir o eixo x com labels especificos
polygon(cord.x,cord.y,col='lightgray') # Comando para sombrear a região no gráficopnorm(1.73)-pnorm(0)pnorm(0)-pnorm(-1.73)1-pnorm(1.73)pnorm(-1.73)pnorm(1.73)-pnorm(0.47)par(mfrow=c(1,3))
regiao=seq(0,1.73,0.01)
cord.x <- c(min(regiao),regiao,max(regiao))
cord.y <- c(0,dnorm(regiao),0)
curve(dnorm(x,0,1),xlim=c(-3,3),main='',xlab="z")
polygon(cord.x,cord.y,col='lightgray')
regiao1=seq(-3,-1.73,0.01)
regiao2=seq(1.73,3,0.01)
cord.x1 <- c(min(regiao1),regiao1,max(regiao1))
cord.x2 <- c(min(regiao2),regiao2,max(regiao2))
cord.y1 <- c(0,dnorm(regiao1),0)
cord.y2 <- c(0,dnorm(regiao2),0)
curve(dnorm(x,0,1),xlim=c(-3,3),main='',xlab="z")
polygon(cord.x1,cord.y1,col='lightgray')
polygon(cord.x2,cord.y2,col='lightgray')
regiao=seq(0.47,1.73,0.01)
cord.x <- c(min(regiao),regiao,max(regiao))
cord.y <- c(0,dnorm(regiao),0)
curve(dnorm(x,0,1),xlim=c(-3,3),main='',xlab="z")
polygon(cord.x,cord.y,col='lightgray')pnorm(5,mean=3,sd=4)-pnorm(2,mean=3,sd=4)par(mfrow=c(1,2))
regiao=seq(2,5,0.01)
cord.x <- c(min(regiao),regiao,max(regiao))
cord.y <- c(0,dnorm(regiao,3,4),0)
curve(dnorm(x,3,4),xlim=c(-9,15),main='X - Normal(3,16)')
polygon(cord.x,cord.y,col='lightgray')
lines(x=c(3,3),y=c(0,dnorm(3,3,4)))
regiao=seq(-0.25,0.5,0.01)
cord.x <- c(min(regiao),regiao,max(regiao))
cord.y <- c(0,dnorm(regiao),0)
curve(dnorm(x,0,1),xlim=c(-3,3),main='Z - Normal Padrão',xlab="z")
polygon(cord.x,cord.y,col='lightgray')
lines(x=c(0,0),y=c(0,fdnorm[x==0]))
regiao=seq(2,5,0.01)
cord.x <- c(min(regiao),regiao,max(regiao))
cord.y <- c(0,dnorm(regiao,3,4),0)
curve(dnorm(x,3,4),xlim=c(-9,15),ylim=c(0,0.4),xlab="",ylab="",xaxs="i",yaxs="i",col="lightblue3",lwd=2)
polygon(cord.x,cord.y,col='lightblue3')
lines(x=c(3,3),y=c(0,dnorm(3,3,4)),col="darkgray",lty=2)
par(new=TRUE)
regiao=seq(-0.25,0.5,0.01)
cord.x <- c(min(regiao),regiao,max(regiao))
cord.y <- c(0,dnorm(regiao),0)
curve(dnorm(x,0,1),xlim=c(-9,15),ylim=c(0,0.4),xlab="",ylab="",xaxs="i",yaxs="i",col="orange2",lwd=2)
polygon(cord.x,cord.y,col='orange2')
lines(x=c(0,0),y=c(0,fdnorm[x==0]),col="gray",lty=2)
legend(5,0.35,c("N(3,16)","N(0,1)"),lwd=c(2,2), col=c("lightblue3","orange2"),pch=c(15,15),bty="n")pnorm(10000,mean=10000,sd=1500)1-pnorm(10000,mean=10000,sd=1500)pnorm(15000,mean=10000,sd=1500)-pnorm(12000,mean=10000,sd=1500)1-pnorm(20000,mean=10000,sd=1500)par(mfrow=c(1,2))
curve(dexp(x,1),xlim=c(0,5),main="f.d.p para X~exponencial(1)",xlab="(a)")
curve(pexp(x,1),xlim=c(0,5),main="f.d.a para X~exponencial(1)",xlab="(b)")Para o cálculo da probabilidade \(P(T>500)\), o comando é:
1-pexp(500,rate=1/500)1-pbinom(6,size=10,prob = 1/2)
1-pnorm(6.5,mean=5, sd=sqrt(2.5))ymax=0.3
barplot(height = dbinom(0:10,size=10,prob = 1/2),col = "white",ylim=c(0,ymax))
par(new=T)
barplot(height = c(rep(0,7),dbinom(7:10,size=10,prob = 1/2)),ylim=c(0,ymax),
border=c(rep(NA,7),rep("black",4)), col = c(rep(NA,7),rep("gray",4)))
par(new=T)
curve(dnorm(x,mean=5, sd=sqrt(2.5)),xlim=c(-0.8,10.8),ylim=c(0.0,ymax),xaxs="i",yaxs="i",ylab="") \(P(3<Y\le 6)=\)
pbinom(6,size=10,prob=1/2)-pbinom(3,size=10,prob=1/2)\(P(3.5\le X \le 6.5)=\)
pnorm(6.5,mean=5,sd=sqrt(2.5))-pnorm(3.5,mean=5,sd=sqrt(2.5))ymax=0.3
barplot(height = dbinom(0:10,size=10,prob = 1/2),col = "white",ylim=c(0,ymax))
par(new=T)
barplot(height = c(rep(0,4),dbinom(4:6,size=10,prob = 1/2),rep(0,4)),ylim=c(0,ymax),
border=c(rep(NA,4),rep("black",3),rep(NA,4)), col = c(rep(NA,4),rep("gray",3),rep(NA,4)))
par(new=T)
curve(dnorm(x,mean=5, sd=sqrt(2.5)),xlim=c(-0.8,10.8),ylim=c(0.0,ymax),xaxs="i",yaxs="i",ylab="") curve(dgamma(x,shape = 3, scale = 1),xlim=c(0,10))par(mfrow=c(1,3))
curve(dchisq(x,df=1),xlim=c(0,10),xlab="(a) df=1")
curve(dchisq(x,df=2),xlim=c(0,10),xlab="(b) df=2")
curve(dchisq(x,df=3),xlim=c(0,10),xlab="(c) df=3")\(Y\sim \chi^2(10)\), então:
\(P(Y>y_0)=0.99\Rightarrow y_0=2.558\).
qchisq(1-.99,df=10)
regiao=seq(2.6,20,0.01)
cord.x <- c(min(regiao),regiao,max(regiao))
cord.y <- c(0,dchisq(regiao,10),0)
curve(dchisq(x,10),xlim=c(0,20),xlab="",ylim=c(0,0.15),ylab="",
col="lightblue3",lwd=2,xaxt="n")
axis(1,at=2.6,labels = "y_0=2.558")
polygon(cord.x,cord.y,col='lightblue3')qchisq(1-0.05,df=10)
regiao=seq(18.3,20,0.01)
cord.x <- c(min(regiao),regiao,max(regiao))
cord.y <- c(0,dchisq(regiao,10),0)
curve(dchisq(x,10),xlim=c(0,20),ylim=c(0,0.15),xlab="",ylab="",
col="lightblue3",lwd=2,xaxt="n")
axis(1,at=18.3,labels = "y_0=18.307")
polygon(cord.x,cord.y,col='lightblue3')
print(paste("Aproximação Normal: ",round(qnorm(1-0.099),3)))
print(paste("Valor /qui-quadrado: ",round(qchisq(1-0.1,df=30),3)))curve(dnorm(x),ylim=c(0,0.5),xlim=c(-3,3),xlab="",col="darkred",ylab="")
par(new=TRUE)
curve(dt(x,df=1),ylim=c(0,0.5),xlim=c(-3,3),xlab="",col="blue",lty=2,ylab="")
legend(-3,0.5,lty=c(1,2), col=c("darkred","blue"),legend=c("normal padrão", "t-student"),bty="n",lwd=c(2,2),cex=0.75)\(P(-1.943 < t(6)< 1.943)=\)
pt(1.943,df = 6)-pt(-1.943,df = 6)\(P( t(6)> 2.447)=\)
1-pt(2.447,df = 6)\(F \sim F(5,7)\), então \(P(F > 3.97)=\)
P_F=pf(3.97,df1 = 5, df2 = 7)
f=qf(0.05,df1 = 5, df2 = 7)print(paste("Q(0.5)=",round(qnorm(0.5),3)))
print(paste("Q(0.25)=",round(qnorm(0.25),3)))
print(paste("Q(0.3)=",round(qnorm(0.3),3)))
print(paste("Q(0.75)=",round(qnorm(0.75),3)))print(paste("Q(0.5)=",round(qexp(0.5,rate = 1/2),3)))pnorm(8.65,mean=10,sd=25)
qnorm(0.8269,mean=10,sd=25)
pexp(0.85,rate=2)
qexp(0.345,rate=2)curve(pnorm(q = x, mean = 0, sd=1), xlim=c(-4,4), col="blue",lwd=2,ylab="",xlab="z")
lines(x=c(0,0), y=c(0,fdanorm[x==0]), lty=2, col="gray")tab<-tabular( estado_civil ~ grau_instrucao*idade_anos*(mean) )
pander(tab)
tab<-tabular( estado_civil*grau_instrucao ~ idade_anos*(mean+sd) )
pander(tab)
tabular( (estado_civil + 1)*(grau_instrucao+1) ~ idade_anos*(mean+sd) )
tab<-tabular( (estado_civil + 1)*(grau_instrucao+1) ~ idade_anos*(mean+sd))
pander(tab)
tabular( (factor(estado_civil) + 1)*(grau_instrucao+1) ~ idade_anos*(mean+sd) )scatter.marginal(tab2_1$idade_anos, tab2_1$salario, ylab="salario", xlab="idade")
scatter.marginal(tab2_1$idade_anos, tab2_1$salario, ylab="salario", xlab="idade",type="boxplot")
scatter.marginal(tab2_1$idade_anos, tab2_1$salario, ylab="salario", xlab="idade",type="dispersion")
dispersao = ggplot(tab2_1, aes(idade_anos,salario)) + geom_point() #Dispersao
ggMarginal(dispersao, type="histogram") # Marginais com histograma
dispersao = ggplot2::ggplot(tab2_1, ggplot2::aes(idade_anos, salario)) + ggplot2::geom_point() #Dispersao
ggMarginal(dispersao,type = "boxplot") # Marginais com boxplot
dispersao = ggplot2::ggplot(tab2_1, ggplot2::aes(idade_anos, salario)) + ggplot2::geom_point() #Dispersao
scatter <- ggplot(tab2_1, aes(idade_anos,salario)) + geom_point() + #dispersão
scale_x_continuous(limits=c(min(idade_anos),max(idade_anos))) + # marginal X
scale_y_continuous(limits=c(min(salario),max(salario))) + # marginal Y
geom_rug(size=0.1) + theme_set(theme_minimal(base_size = 18))
scatter
ggMarginal(dispersao,type = "density") # Marginais com função densidade empíricapairs( ~ n_filhos + salario + idade_anos, pch=16, col="blue")
pairs( ~ n_filhos + salario + idade_anos, pch=16, col="blue",upper.panel = NULL)contour(x=1:nrow(volcano), y=1:ncol(volcano), volcano, col="darkred",xlab = "Sul - Norte", ylab = "Leste - Oeste")
title("Mapa topográfico: Maunga Whau", font = 4) # Titulo do GráficoFigura 8.A: Gráfico de curvas de nível do vulcão Maunga Whau em função da latitude e longitude.
library(colorRamps) # Paleta de cores que contem a sequencia de cores `blue2red`.
filled.contour(x=1:nrow(volcano), y=1:ncol(volcano), z=volcano, color = blue2red,
plot.title = title(main = "Mapa topográfico: Maunga Whau \n ",
xlab = "Sul - Norte", ylab = "Leste - Oeste"),
plot.axes = { axis(1, seq(100, 800, by = 100))
axis(2, seq(100, 600, by = 100)) },
key.axes = axis(4, seq(90, 190, by = 10)))# outra sugestão de cor: `color = terrain.colors` ou não declarar esta opção para visualizar a cor padrão. cov(salario,idade_anos)
cov(idade_anos,salario)Sigma=cbind(c(1,0.6),c(0.6,1)) # matriz de variancias e covariancias
mu=c(0,0) # Vetor de médias
Z.sim<-rmvnorm(n = 1000, mean=mu, sigma=Sigma) # simula 100 observações de uma normal bivariada
fdamvnorm<-function(upper,lower=c(-Inf,-Inf),mu=rep(0,times=length(upper)),Sigma=diag(length(upper))){
return(pmvnorm(lower=c(-Inf,-Inf),upper=upper,mean=mu,sigma=Sigma)[[1]]) # Calcula a integral dupla da fdp da normal multivariada entre no retangulo `lower X upper`.
}
fdamvnorm(upper=c(0,0)) # FDA considerando o terceiro quadrante
fdamvnorm(upper=c(Inf,0)) # FDA considerando o terceiro e quarto quadrantes
fdamvnorm(c(0,Inf)) # FDA considerando o segundo e terceiro quadrantes
fdamvnorm(lower=c(-1.96,-1.96),upper=c(1.96,1.96)) #FDA considerando o quadrado (-1.96,-1.96)X(1.96,1.96)
Sigma=cbind(c(1,0.6),c(0.6,1)) # matriz de variancias e covariancias
mu=c(0,0) # Vetor de médias
dmvnorm(x=c(0,0), mean=mu, sigma=Sigma)xy<-seq(-4,4,0.1)
mu=c(0,0)
Sigma=diag(2)
Z1<-matrix(NA,length(xy),length(xy))
for(i in 1:length(xy)){
for(j in 1:length(xy)){
Z1[i,j]<-dmvnorm(x=c(xy[i],xy[j]), mean=mu, sigma=Sigma)
}
}persp(Z1, phi = 30, theta = -45,col=heat.colors(ncol(Z1)*nrow(Z1),.8), ticktype = "detailed",xlab="X")
# `phi` e `theta` são os parâmetros de deslocamento da perspectiva da figura.
# heat.colors(.) modifica a paleta de cores do gráfico
# box=FALSE elimina a caixa que envolve o gráfico
# scale =FALSE define que as variáveis não devem ser padronizadas
# ticktype: tipo de eixos. Pode assumir valores "simple" ou "detailed"xy<-seq(-4,4,0.1)
mu=c(0,0)
Sigma=cbind(c(1,0.6),c(0.6,1))
Z2<-matrix(NA,length(xy),length(xy))
for(i in 1:length(xy)){
for(j in 1:length(xy)){
Z2[i,j]<-dmvnorm(x=c(xy[i],xy[j]), mean=mu, sigma=Sigma)
}
}
persp(Z2, phi = 30, theta = -45,col=heat.colors(ncol(Z2)*nrow(Z2),.8), ticktype = "detailed", xlab="X")par(mfrow=c(1,2))
contour(x=1:nrow(Z2), y=1:ncol(Z2), Z2, col="darkred",xlab = "X", ylab = "Y")
contour(x=1:nrow(Z1), y=1:ncol(Z1), Z1, col="darkred",xlab = "X", ylab = "Y")u<-runif(10,0,1) # gera 10 observações aleatórias de uma v.a. uniforme[0,1]
u
u <- runif(1000,0,1) # gera 1000 observações aleatórias de uma v.a. uniforme[0,1]
x <- qnorm(u,mean=0, sd = 1) # Calcula os quantis para o vetor u simulado da uniforme.par(mfrow=c(1,2))
hist(u, freq=FALSE, main="Histograma da amostra da \n distribuição Uniforme simulada", col="lightblue4",border="gray")
hist(x, freq=FALSE, main="Histograma da variável X simulada a partir \n do resultado do Teorema 9.1", col="lightblue4",border="gray")u <- runif(1000,0,1) # gera 1000 observações aleatórias de uma v.a. uniforme[0,1]
x <- sqrt(u) # Calcula a f.d.a. inversa para o vetor u gerado.
par(mfrow=c(1,2))
hist(u, freq=FALSE, main="Histograma da amostra da \n distribuição Uniforme simulada", col="lightblue4",border="gray")
hist(x, freq=FALSE, main="Histograma da variável x simulada a partir \n do resultado do Teorema 9.1", col="lightblue4",border="gray")par(mfrow=c(1,2))
x<-seq(0,1.75,0.01) # Produzindo o vetor da abscissa
f_x<-c(2*x[x<=1],0*x[x>1]) # Calculando f(x)
plot(x,f_x, col='darkred', ylab="f(x)",xaxt="n",type="l") #Grafico de f(x)
lines(x=c(0,1),y=c(2,2), col="gray",lty=2)
F_x<-c((x[x<=1])^2,rep(1,length(x))[x>1]) # Calculando F(x)
plot(x,F_x, col='darkblue', ylab="F(x)",xaxt="n",type="l") #Grafico de F(x)
lines(c(0,1),c(1,1), col="gray",lty=2) # Adicionando linas explicativas tracejadas
lines(c(1,1),c(0,1), col="gray",lty=2)
lines(c(0,.71),c(0.51,.51), col="gray",lty=2)
lines(c(.71,.71),c(0,.51), col="gray",lty=2)
text(0,.51,"u=0.5", pos=4) # inserindo texto no grafico
text(0.71,0,"x=0.71")rbernoulli(p=0.52)
u<-rbernoulli(n=10,p=0.52)
u
aa<-rbernoulli(n=1000, p=0.52)
sum(aa)/1000xbinomial<-numeric()
for( i in 1:20){
xbinomial[i]<-sum(rbernoulli(n=10,p=0.52))
}
xbinomial
rbinom(n=20, size = 10, p=0.52)u<-runif(n=5,0,1)
qexp(u,rate=1/2)
rexp(n=5,rate=1/2)qnorm(0.230,mean=10, sd = .4) # Calcula os quantis para a normal com média 10 e sd = 0.4x<-rbinom(20,10,0.5)
x
y<-rpois(20,1.7)
yx<-rbinom(20,10,0.5)
y<-rpois(20,1.7)
z<-runif(100,0,1)
b<-rbinom(15,1,0.7)par(mfrow=c(2,2))
hist(x, col="lightblue4",border="white")
hist(y, col="lightblue4",border="white")
hist(z, col="lightblue4",border="white")
hist(b, col="lightblue4",border="white")z<-rnorm(500,0,1)
y<-rnorm(200,10,3)
t<-rt(500,35)
Exp<-rexp(500,2)
w<-rchisq(300,5)
f<-rf(500,10,12)``{r, eval=FALSE} par(mfrow=c(3,2)) hist(z, col=“lightblue4”,border=“white”) hist(y, col=“lightblue4”,border=“white”) hist(t, col=“lightblue4”,border=“white”) hist(Exp, col=“lightblue4”,border=“white”) hist(w, col=“lightblue4”,border=“white”) hist(f, col=“lightblue4”,border=“white”) ```
u1<-runif(200,0,1)
u2<-runif(200,0,1)plot(u1,u2, pch=20, col="darkblue") u1<-runif(10000,0,1)
u2<-runif(10000,0,1)
P<-data.frame(u1,u2)
x<-qnorm(u1)
y<-qnorm(u2)scatter.marginal(x,y,breaks=30) # Função criada no Capitulo 8 para construir este tipo de gráficolibrary(mvtnorm)
library(scatterplot3d)
scatterplot3d(x,y, dmvnorm(x=cbind(x,y)), highlight.3d=TRUE, pch=16, zlab="Z")contour(interp(x,y, dmvnorm(x=cbind(x,y)))) ##### Amostragem sem reposição: Utilize a opção 'replace = FALSE'
AASs = sample(x = 1:dim(tab2_1)[1], size=15, replace = FALSE )
# Elementos selecionados na amostra
AASs
##### Amostragem com reposição: Utilize a opção 'replace = TRUE'
AASc = sample(x = 1:dim(tab2_1)[1], size=15, replace = TRUE )
# Elementos selecionados na amostra
AASc
print(kable(tab2_1[AASs,], caption = '**Tabela 10.A**: Amostra Aleatória Simples sem Reposição da Tabela 2.1'))
print(kable(tab2_1[AASc,], caption = '**Tabela 10.B**: Amostra Aleatória Simples com Reposição da Tabela 2.1'))x<-rnorm(5,mean=167, sd= 5)
xx_normal<-rnorm(10000,mean=167, sd= 5)
hist(x_normal, col="lightblue4", border="white",freq = FALSE)
#Inicializando as variaveis como vetores numericos
## Media e variancia para amostras de tam 15
xbar15<-numeric()
var_amostral15<-numeric()
## Media e variancia para amostras de tam 30
xbar30<-numeric()
var_amostral30<-numeric()
## Media e variancia para amostras de tam 100
xbar100<-numeric()
var_amostral100<-numeric()
for ( i in 1:200){
# Extraindo amostras de tamanho 15 e calculando a média e variancia
smp<-sample(x_normal,size = 15)
xbar15[i]<-mean(x_normal[smp])
var_amostral15[i]<-var(x_normal[smp])
# Extraindo amostras de tamanho 30 e calculando a média e variancia
smp<-sample(x_normal,size = 30)
xbar30[i]<-mean(x_normal[smp])
var_amostral30[i]<-var(x_normal[smp])
# Extraindo amostras de tamanho 100 e calculando a média e variancia
smp<-sample(x_normal,size = 100)
xbar100[i]<-mean(x_normal[smp])
var_amostral100[i]<-var(x_normal[smp])
}par(mfrow=c(2,3))
hist(xbar15, col="lightblue4", border="white",freq = FALSE, breaks = 15, main="Média com n=15")
hist(xbar30, col="lightblue4", border="white",freq = FALSE, breaks = 15, main="Média com n=30")
hist(xbar100, col="lightblue4", border="white",freq = FALSE, breaks = 15, main="Média com n=100")
hist(var_amostral15, col="lightblue4", border="white",freq = FALSE, breaks = 15, main="Variância com n=15")
hist(var_amostral30, col="lightblue4", border="white",freq = FALSE, breaks = 15, main="Variância com n=30")
hist(var_amostral100, col="lightblue4", border="white",freq = FALSE, breaks = 15, main="Variância com n=100")plot (density(x_normal), ylim=c(0,1), xlab="", main = "Funções densidades empíricas")
lines (density(xbar15),col="red")
lines (density(xbar30),col="blue")
lines (density(xbar100),col="orange")
legend("topright", c("X~N(167,25)", "Média com n=15", "Média com n=30", "Média com n=100"),
col=c("black", "red", "blue", "orange"), lty=c(1,1,1,1))## Mediana para amostras de tam 5
md5<-numeric()
for ( i in 1:200){
# Extraindo amostras de tamanho 5 e calculando a mediana
smp<-sample(x_normal,size = 5)
md5[i]<-median(x_normal[smp])
}
print(summary(md5))
cat("E(md)=",mean(md5), ", Var(md)=",var(md5),", dp(md)=",sd(md5),", min(md)=",min(md5),", max(md)=",max(md5))hist(md5, col="lightblue3", border="white",freq = FALSE, main="Mediana com n=5")
# breaks=seq(min(md5)-1,max(md5)+1,2),n=cbind(c(2, 5,25), # numero de amostras para os gráficos da primeira linha
c(2,10,30), # numero de amostras para os gráficos da segunda linha
c(3, 5,10), # numero de amostras para os gráficos da terceira linha
c(5,15,25)) # numero de amostras para os gráficos da quarta linha
X<-list(
x1=matrix(NA,nrow=200,ncol=3), # Matriz que armazenara as médias da primeira linha
x2=matrix(NA,nrow=200,ncol=3), # Matriz que armazenara as médias da segunda linha
x3=matrix(NA,nrow=200,ncol=3), # Matriz que armazenara as médias da terceira linha
x4=matrix(NA,nrow=200,ncol=3)) # Matriz que armazenara as médias da quarta linha
for ( i in 1:3){ # iterando ao longo do número de amostras de cada linha
for ( j in 1: 200){ # Construindo 200 observacoes para cada gráfico
X$x1[j,i] = mean( runif( n=n[i,1])) # Calculando a média pra cada uma das 200 simulacoes de uma distribuicao uniforme
X$x2[j,i] = mean( rexp( n=n[i,2], rate=2)) # Calculando a média pra cada uma das 200 simulacoes de uma distribuicao exponencial
X$x3[j,i] = mean(rchisq( n=n[i,3], df=1)) # Calculando a média pra cada uma das 200 simulacoes de uma distribuicao qui-quadrado
X$x4[j,i] = mean( rbeta( n=n[i,4], shape1 = .2, shape2 = .2)) # Calculando a média pra cada uma das 200 simulacoes de uma distribuicao beta
}
}
par(mfrow=c(4,4))
plot(dunif(x = seq(0,1,0.05)),type="l", ylab="f(x)", main="fdp-Uniforme(0,1)", xlab="X")
hist(X$x1[,1], col = "lightblue3", border = "white", main = paste("n = ", n[1,1]), freq = FALSE, xlab="X_bar")
hist(X$x1[,2], col = "lightblue3", border = "white", main = paste("n = ", n[2,1]), freq = FALSE, xlab="X_bar")
hist(X$x1[,3], col = "lightblue3", border = "white", main = paste("n = ", n[3,1]), freq = FALSE, xlab="X_bar")
plot(dexp(x = seq(0,1,0.05),rate=2),type="l", ylab="f(x)", main="fdp-Exponencial(1/2)", xlab="X")
hist(X$x2[,1], col = "lightblue3", border = "white", main = paste("n = ", n[1,2]), freq = FALSE, xlab="X_bar")
hist(X$x2[,2], col = "lightblue3", border = "white", main = paste("n = ", n[2,2]), freq = FALSE, xlab="X_bar")
hist(X$x2[,3], col = "lightblue3", border = "white", main = paste("n = ", n[3,2]), freq = FALSE, xlab="X_bar")
plot(dchisq(x = seq(0,1,0.05),df=1),type="l", ylab="f(x)", main="fdp-QuiQuad(1)", xlab="X")
hist(X$x3[,1], col = "lightblue3", border = "white", main = paste("n = ", n[1,3]), freq = FALSE, xlab="X_bar")
hist(X$x3[,2], col = "lightblue3", border = "white", main = paste("n = ", n[2,3]), freq = FALSE, xlab="X_bar")
hist(X$x3[,3], col = "lightblue3", border = "white", main = paste("n = ", n[3,3]), freq = FALSE, xlab="X_bar")
plot(dbeta(x = seq(0,1,0.05), shape1 = .2, shape2 = .2),type="l", ylab="f(x)", main="fdp-Beta(0.2, 0.2)", xlab="X")
hist(X$x4[,1], col = "lightblue3", border = "white", main = paste("n = ", n[1,4]), freq = FALSE, xlab="X_bar")
hist(X$x4[,2], col = "lightblue3", border = "white", main = paste("n = ", n[2,4]), freq = FALSE, xlab="X_bar")
hist(X$x4[,3], col = "lightblue3", border = "white", main = paste("n = ", n[3,4]), freq = FALSE, xlab="X_bar")x_barr<-sapply(1:200, function(i){return(mean(rnorm(n=100,mean=500,sd=10)))})
hist(x_barr, col = "lightblue3", border = "white", freq=FALSE)
sum( x_barr >= 498 & x_barr <= 502 )/200x_barr <- numeric()
x_var <- numeric()
for (i in 1:200){
smp <- rbinom(n = 10, size = 1, prob = 0.3)
x_barr[i] <- mean(smp)
x_var[i] <- var(smp)
}
p_hat = sum( x_barr ) / 200
p_hat
hist(x_barr, col="lightblue3", border="white", breaks=seq(0,1,0.02))
x_barr <- numeric()
x_var <- numeric()
for (i in 1:200){
smp <- rbinom(n = 100, size = 1, prob = 0.3)
x_barr[i] <- mean(smp)
x_var[i] <- var(smp)
}
p_hat = sum( x_barr ) / 200
p_hat
hist(x_barr, col="lightblue3", border="white", breaks=seq(0,1,0.02))x<-c(1,3,5,5,7)
smp<-matrix(NA,nrow=200,ncol=2)
for ( i in 1:200){
smp[i,] <- sample(x, size=2, replace=TRUE)
}
hist(smp, col="lightblue3", border="white")
x<-c(1,3,5,5,7)
S2<-numeric()
md<-numeric()
var_hat<-numeric()
for ( i in 1:200){
smp<-sample(x, size=3, replace=TRUE)
S2[i] <- var(smp)
var_hat[i] <- var(smp)*2/3
md[i] <- median(smp)
}hist(S2, col="lightblue3", border="white", main="Figura 10.6")
hist(md, col="lightblue3", border="white", main="Figura 10.7")
hist(var_hat, col="lightblue3", border="white", main="Figura 10.8")x<-c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)
sample(x, 7)
sample(x, 7, replace=T)tab11.1<- data.frame(
x=c(1.2,1.5,1.7,2,2.6),
y=c(3.9,4.7,5.6,5.8,7.0))
mean(tab11.1$x)
mean(tab11.1$y)
plot(tab11.1$x,tab11.1$y,pch=16,col= "darkblue")
tab11.1$x_3=3*tab11.1$x
tab11.1$y_3x=tab11.1$y - tab11.1$x_3
tab11.1$y_3x_2=tab11.1$y_3x^2
print(kable(tab11.1,col.names = c("$X$","$Y$", "$3X$","$Y-3X$","$(Y-3X)^2$"), caption = "**Tabela 11.1**: Análise do modelo $Ŷ=3X$"))
print(paste("Total da soma de quadrados dos resíduos(sum(tab11.1$y_3x_2)):",sum(tab11.1$y_3x_2)))
theta_hat=sum(tab11.1$x*tab11.1$y)/sum(tab11.1$x^2)x_sim <- rexp(1000, rate = 2) # 1000 observações amostradas de um modelo exponencial de média 1/2
mean(x_sim)
hist(x_sim, col="lightblue3", border="white", freq = FALSE)
MV <- function(rate) {
R = dnorm(x_sim, rate)
return(-sum(log(R))) # devolvendo '-log(verossimilhanca)'
}
emv_rate_x <- mle(MV, # MV é a função que se deseja otimizar
start = list(rate=1)) # start é uma lista com 'chutes' iniciais para o processo de otimização.
emv_rate_xsamples=20 # Numero de amostras
N=25 # Numero de observações por amostra
x_sim<-matrix(NA,N,samples) # x_sim armazenara as amostras simuladas
for( i in 1:samples){
x_sim[,i]<-rnorm(n = N, mean= 5, sd=3) # Simulando 20 N(5,3^2)
}
cont=0
for( i in 1:samples){
llim=mean(x_sim[,i])-1.96*sd(x_sim[,i])/sqrt(N) # limite inferior do IC da i-ésima amostra
ulim=mean(x_sim[,i])+1.96*sd(x_sim[,i])/sqrt(N) # limite inferior do IC da i-ésima amostra
if (llim>5 || ulim<5){
line_col="red"
cont=cont+1
}
else line_col="black"
if(i!=1)
par(new=T)
plot(x=i,y=mean(x_sim[,i]), xlim=c(0,21), ylim=c(0,10), pch=20,cex=1, col=line_col, ylab="Médias", xlab="Amostras") # plotando as médias de cada amostra
lines(x = c(i,i), y =c(llim,ulim), lwd=1, col=line_col) # Adicionando as linhas dos ICs de cada amostra
abline(h=5,lty = 'dashed', col="blue") # Valor real da média
}IC(n=25, media = 485, desvpad=sqrt(100), gama=0.95)IC(n=400, media = 0.6, desvpad=sqrt(1/4), gama=0.95)p=80/400
dp=sqrt(p*(1-p))
IC(n=400, media = p, desvpad=dp, gama=0.9)
# IC Conservador
p=1/2
dp=sqrt(p*(1-p))
IC(n=400, media = 80/400, desvpad=dp, gama=0.9)amostra_original<-rexp(50, rate=3) # Vamos simular uma exponencial de média 1/3
n=length(amostra_original)
md=numeric() # vetor que armazenará as medianas calculadas ao longo das amostras bootstrap
b=30 # Número de amostras bootstrap a serem resorteadas baseadas na amostra original.
for( i in 1:b){ # Vamos calcular 30 valores da mediana para construir sua dist. empirica
boot<-sample(amostra_original, size = n, replace = TRUE) # Tomando amostras de tamanho n, com reposicao
md[i]<-median(boot)
}
hist(md, col="lightblue3", border= "white")
ep_md=numeric() # vetor que armazenará os erros padrões da mediana para cada rodada do bootstrap
md=numeric() # vetor que armazenará as medianas calculadas ao longo das amostras bootstrap
b=c(30,50,100,500,1000) # Número de amostras bootstrap a serem resorteadas baseadas na amostra original.
for (j in 1:5){ # neste for, percorreremos os valores de b
for( i in 1:b[j]){ # Vamos calcular 30 valores da mediana para construir sua dist. empirica
boot<-sample(amostra_original, size = n, replace = TRUE) # Tomando amostras de tamanho n, com reposicao
md[i]<-median(boot)
}
ep_md[j]=sd(md)/b[j]
}
plot(1:5, ep_md, cex=2, type="b", ylab="Erro Padrão", xlab="b: tamanho do bootstrap")amostra_original<-c(2,5,3,4,6)
n=length(amostra_original)
b=5
boot<-matrix(NA,n,b) # armazena as amostras geradas
md=numeric()
xbarr=numeric()
for( i in 1:b){
boot[,i]<-sample(amostra_original, size = n, replace = TRUE)
md[i]<-median(boot[,i])
xbarr[i]<-mean(boot[,i])
}tab11.2<-data.frame(x1_boot=boot[,1],
x2_boot=boot[,2],
x3_boot=boot[,3],
x4_boot=boot[,4],
x5_boot=boot[,5],
md=md,
media=xbarr)
print(kable(tab11.2,col.names = c("Bootstrap:x1*","Bootstrap:x2*","Bootstrap:x3*","Bootstrap:x4*","Bootstrap:x5*", "md(x*)","média(x*)"),caption = "**Tabela 11.2**: Procedimento *bootstrap*."))
### Calculando agora o EP para a mediana a partir das amostras geradas:
ep=0
for( i in 1:b){
ep=ep+(md[i]-mean(md))^2/(b-1)
}
ep=sqrt(ep)
ep### Calculando agora o EP para a mediana a partir das amostras geradas:
ep_media=0
for( i in 1:b){
ep_media=ep_media+(xbarr[i]-mean(xbarr))^2/(b-1)
}
ep_media=sqrt(ep_media)
ep_media
IC(n = 1,media = 5.2, desvpad = ep, gama = 0.95)
amostra_original<-c(2,5,3,4,6)
n=length(amostra_original)
b=50
boot<-matrix(NA,n,b) # armazena as amostras geradas
md=numeric()
xbarr=numeric()
for( i in 1:b){
boot[,i]<-sample(amostra_original, size = n, replace = TRUE)
md[i]<-median(boot[,i])
xbarr[i]<-mean(boot[,i])
}
### Calculando o EP para a mediana a partir das 50 amostras geradas:
ep=0
for( i in 1:b){
ep=ep+(md[i]-mean(md))^2/(b-1)
}
ep=sqrt(ep)
ep
### Calculando o EP para a media a partir das50 amostras geradas:
ep_media=0
for( i in 1:b){
ep_media=ep_media+(xbarr[i]-mean(xbarr))^2/(b-1)
}
ep_media=sqrt(ep_media)
ep_mediabootstrapamostra_original<-c(2,5,3,4,6)
### Exemplo 11.19: Mediana
boot_11_19<-bootstrap(x=amostra_original, nboot=5, theta=function(x){median(x)})
boot_11_19$thetastar
### Exemplo 11.19: Mediana
boot_11_20<-bootstrap(x=amostra_original, nboot=5, theta=function(x){mean(x)})
boot_11_20$thetastar
#Se desejarmos, como nos exemplos, a estimativa do erro padrão, faríamos:
err_pad<- function(x){sqrt(var(x)/length(x))}
bootstrap(x=amostra_original, nboot=15, theta=function(x){median(x)}, func=err_pad)
bootstrap(x=amostra_original, nboot=15, theta=function(x){mean(x)}, func=err_pad)qnorm(p = 0.05, mean=155, sd =4)x<-seq(125,167,0.1) # eixo X do gráfico
regiao1=seq(148.42,167,0.01) # Regiao a ser sombreada na figura
cord.x <- c(min(regiao1),regiao1,max(regiao1)) # coordenadas X da regiao
cord.y <- c(0,dnorm(regiao1,mean=145,sd=5.76),0) # coordenadas Y da regiao
curve(dnorm(x,mean=145,sd=5.76),xlim=c(125,167),ylim=c(0,0.12),xlab="",ylab="",xaxs="i",yaxs="i",col="black",lwd=2, xaxt='n') # Curva da distribuicao normal
polygon(cord.x,cord.y,col='orange2') # plotando a regiao
lines(x=c(145,145),y=c(0,dnorm(145,145,5.76)),col="gray",lty=2) # linha vertical da media
par(new=TRUE) # inserindo a segunda curva
regiao2=seq(125,148.42,0.01) # Regiao a ser sombreada na figura
cord.x <- c(min(regiao2),regiao2,max(regiao2))
cord.y <- c(0,dnorm(regiao2,mean=155,sd=4),0)
curve(dnorm(x,mean=155,sd=4),xlim=c(125,167),ylim=c(0,0.12),xlab="",ylab="",xaxs="i",yaxs="i",col="black",lwd=2, xaxt='n')
polygon(cord.x,cord.y,col='lightblue3')
lines(x=c(155,155),y=c(0,dnorm(155,155,4)),col="darkgray",lty=2)
axis(side=1,at = c(145,148.42,155), labels = c("145","148.42","155")) # 'labels' do eixo X
legend("topleft", legend=c("RC=5%", "Poder=7.93%"), pch=c(15,15), col=c("lightblue3","orange2"), cex=2, bty = "n") # legendax<-seq(125,167,0.1)
regiao1=seq(125,148.42,0.01) # Regiao a ser sombreada na figura
cord.x <- c(min(regiao1),regiao1,max(regiao1))
cord.y <- c(0,dnorm(regiao1,mean=155,sd=4),0)
curve(dnorm(x,mean=155,sd=4),xlim=c(125,167),ylim=c(0,0.15),
xlab="",ylab="",xaxs="i",yaxs="i",col="black",lwd=2, xaxt='n')
polygon(cord.x,cord.y,col='lightblue3')
lines(x=c(155,155),y=c(0,dnorm(155,155,4)),col="darkgray",lty=2)
axis(side=1,at = c(148.42,155), labels = c("148.42","mu=155"))
par(new=TRUE)
# Curva tracejada
curve(dnorm(x,mean=148.42,sd=3),xlim=c(125,167),ylim=c(0,0.15),
xlab="",ylab="",xaxs="i",yaxs="i",col="gray",lty=2, lwd=2, xaxt='n')
par(new=TRUE)
# Curva tracejada
curve(dnorm(x,mean=145,sd=3),xlim=c(125,167),ylim=c(0,0.15),
xlab="",ylab="",xaxs="i",yaxs="i",col="gray",lty=2, lwd=2, xaxt='n') x<-seq(130,190,0.1)
regiao2=seq(130,147.16,0.01)
regiao1=seq(162.84,190,0.01)
cord.x1 <- c(min(regiao1),regiao1,max(regiao1))
cord.x2 <- c(min(regiao2),regiao2,max(regiao2))
cord.y1 <- c(0,dnorm(regiao1,mean=155,sd=4),0)
cord.y2 <- c(0,dnorm(regiao2,mean=155,sd=4),0)
curve(dnorm(x,mean=155,sd=4),xlim=c(125,190),ylim=c(0,0.15),
xlab="",ylab="",xaxs="i",yaxs="i",lwd=2, xaxt='n')
polygon(cord.x1,cord.y1,col='lightblue3')
polygon(cord.x2,cord.y2,col='lightblue3')
lines(x=c(155,155),y=c(0,dnorm(155,155,4)),col="darkgray",lty=2)
axis(side=1,at = c(147.16,155,162.84),
labels = c("X_c1","155","X_c2"))
par(new=TRUE)
# Curva tracejada
curve(dnorm(x,mean=145,sd=4.5),xlim=c(130,190),ylim=c(0,0.15),
xlab="",ylab="",xaxs="i",yaxs="i",col="gray",lty=2, lwd=2, xaxt='n')
par(new=TRUE)
# Curva tracejada
curve(dnorm(x,mean=140,sd=3),xlim=c(130,190),ylim=c(0,0.15),
xlab="",ylab="",xaxs="i",yaxs="i",col="gray",lty=2, lwd=2, xaxt='n')
par(new=TRUE)
# Curva tracejada
curve(dnorm(x,mean=170,sd=3),xlim=c(130,190),ylim=c(0,0.15),
xlab="",ylab="",xaxs="i",yaxs="i",col="gray",lty=2, lwd=2, xaxt='n')
par(new=TRUE)
# Curva tracejada
curve(dnorm(x,mean=175,sd=4.5),xlim=c(130,190),ylim=c(0,0.15),
xlab="",ylab="",xaxs="i",yaxs="i",col="gray",lty=2, lwd=2, xaxt='n') RC(tipo = "bilateral",480,520,500,5,0.01)pc<-RC(tipo = "inferior", x_min = 0.45,x_max = 0.75, media = 0.60, desvpad = sqrt(0.24/200), alpha = 0.05)limites<-RC(tipo="bilateral", alpha=0.01, media=500, desvpad=5, grafico = FALSE)mu_abaixo=c(500,498,495,492,490,487,485,480,475)
mu_acima=c(500,502,505,508,510,513,515,520,525)
tab12.1<-data.frame(mu_abaixo,mu_acima,
pi_mu=(1-sapply(mu_acima,function(mu){return(pnorm(limites[2],mean=mu,sd=5)-pnorm(limites[1],mean=mu,sd=5))}))*100,
beta_mu=sapply(mu_acima,function(mu){return(pnorm(limites[2],mean=mu,sd=5)-pnorm(limites[1],mean=mu,sd=5))})*100)
print(kable(tab12.1,col.names = c("À esquerda de 500", "À direita de 500", "pi(mu) (em %)","beta(mu) (em %)"),digits = 1))par(mfrow=c(6,1))
RC_(tipo="bilateral", alpha=0.01, media=500, desvpad=5, x_min=460,x_max=540,media_=500)
RC_(tipo="superior" , alpha=0.005, media=502, desvpad=5, x_min=460,x_max=540,media_=500)
RC_(tipo="superior" , alpha=0.005, media=510, desvpad=5, x_min=460,x_max=540,media_=500)
RC_(tipo="superior" , alpha=0.005, media=515, desvpad=5, x_min=460,x_max=540,media_=500)
RC_(tipo="inferior" , alpha=0.005, media=495, desvpad=5, x_min=460,x_max=540,media_=500)
RC_(tipo="inferior" , alpha=0.005, media=487, desvpad=5, x_min=460,x_max=540,media_=500) x=c(475,480,485,487,490,492,495,498,500,502,505,508,510,513,515,520,525)
plot(x,(1-sapply(x,function(mu){return(pnorm(limites[2],mean=mu,sd=sqrt(400/16))-pnorm(limites[1],mean=mu,sd=sqrt(400/16)))}))*100,type="b", ylab="", xlim=c(475,525),ylim=c(0,100))
par(new=T)
limites_100<-RC(tipo="bilateral", alpha=0.01, media=500, desvpad=sqrt(400/100), grafico = FALSE)
plot(x,(1-sapply(x,function(mu){return(pnorm(limites_100[2],mean=mu,sd=sqrt(400/100))-pnorm(limites_100[1],mean=mu,sd=sqrt(400/100)))}))*100,type="l", ylab="", lty=2, xlim=c(475,525), ylim=c(0,100))limites_100<-RC(tipo="bilateral", alpha=0.01, media=500, desvpad=sqrt(400/100), grafico = FALSE)val_p=pnorm(0.52, 0.6,sqrt(0.24/200))RC(tipo = "inferior", x_min = 0.5,x_max=0.7,media = 0.6,desvpad = sqrt(0.24/200),alpha = val_p)(1-pnorm(314,300,30/sqrt(10)))*2#Calculo da estatistica observada:
q_obs=15*169/100
q_obs
#Calculo da RC:
alpha=0.05 # nivel de significancia desejado
n=16 # numero de pacotes
qc_inf=qchisq(p = alpha/2,df=n-1) #valor critico da regiao inferior
qc_sup=qchisq(p = 1- alpha/2,df=n-1) #valor critico da regiao superior
#Valores críticos:
c(qc_inf,qc_sup)# Plotando o gráfico:
x_min=0
x_max=35
x<-seq(x_min,x_max,0.1)
regiao1=seq(0,qchisq(p=alpha/2, df=n-1),0.001)
regiao2=seq(qchisq(p=1-alpha/2, df=n-1),x_max,0.001)
cord.x1 <- c(min(regiao1),regiao1,max(regiao1))
cord.x2 <- c(min(regiao2),regiao2,max(regiao2))
cord.y1 <- c(0,dchisq(regiao1,df=n-1),0)
cord.y2 <- c(0,dchisq(regiao2,df=n-1),0)
curve(dchisq(x,df=n-1),xlim=c(x_min,x_max),xlab="",ylab="",xaxs="i",yaxs="i",lwd=2, xaxt='n')
polygon(cord.x1,cord.y1,col='orange2')
polygon(cord.x2,cord.y2,col='orange2')
axis(side=1,at = c(qchisq(p=alpha/2,df=n-1), qchisq(p=1-alpha/2,df=n-1)),
labels = round(c(qchisq(p=alpha/2,df=n-1), qchisq(p=1-alpha/2,df=n-1)),3))x<-c(253,187,96,450,320,105)
s_o<-var(x)
#Variancia observada:
s_o
alpha=0.1 # nivel de significancia desejado
n=6
qc_inf=qchisq(p = alpha/2,df=n-1) #valor critico da regiao inferior
qc_sup=qchisq(p = 1- alpha/2,df=n-1) #valor critico da regiao superior
#Valores críticos:
c(qc_inf,qc_sup)
pchisq(qc_sup,df=n-1)-pchisq(qc_inf,df=n-1)t_c=qt(p=1-0.05,df=24)
t_c
# Valor observado da estatística t:
T=sqrt(25)*(31.5-30)/3
T
1-pt(T,df=24)
# IC para mu
alpha=0.05
n=25
lim_inf=31.5+qt(alpha/2,df=n-1)*(3/sqrt(n))
lim_sup=31.5+qt(1-alpha/2,df=n-1)*(3/sqrt(n))
print(paste("IC(mu;",1-alpha,")=[",round(lim_inf,2),";",round(lim_sup,2),"]"))maqA<-c(145,127,136,142,141,137)
maqB<-c(143,128,132,138,142,132)
s2_a<-var(maqA)
s2_b<-var(maqB)
teste_13_2<-var.test(x = maqA, y = maqB, conf.level = 0.9, alternative = "two.sided")
teste_13_2teste_var(s1=85,s2=8,n1 = 6,n2=6,alpha = 0.1)teste_t(x_barra1=68,x_barra2=76,s1=50,s2=52,n1=12,n2=15,alpha=0.05, H_1="A<B",var_iguais=TRUE)
IC_diff_medias(x_barra1=68,x_barra2=76,s1=50,s2=52,n1=12,n2=15,gama=0.95)teste_t(x_barra1=70.5,x_barra2=84.3,s1=81.6,s2=210.8,
n1=15,n2=20,alpha=0.05, H_1="A!=B", var_iguais=FALSE)pwilcox(q=87-10*11/2,m=10,n=10)controle<-c(1.3, 1.5, 2.1)
tratamento<-c(1.5,2.5)
pwilcox(q=7.5,m=3,n=2)
wilcox.test(controle,tratamento, alternative = "less")T<-c(0.6,0.63,0.83,0.85,0.91,0.95,1.01,1.03,1.03,1.16,
1.19,1.2,1.26,1.28,1.3,1.37,1.45,1.54,1.68,2.2)
C<-c(0.52,0.77,0.79,0.79,0.81,0.81,0.89,0.98,1.01,1.18,
1.19,1.2,1.34,1.36,1.38,1.43,1.64,1.71,2.16,2.25)
hist(T,breaks=seq(0,2.8,0.4),col="lightgray", main="", freq = FALSE,yaxt='n',ylab="")hist(C,breaks=seq(0,2.8,0.4),col="lightgray", main="", freq = FALSE,yaxt='n',ylab="")media<-(c(C,T))
tipo<-factor(c(rep("C",20),rep("T",20)))
posto<-rank(c(C,T))
tab13_7<-data.frame(media,tipo,posto)
tab13_7<-tab13_7[order(tab13_7$posto),]
print(kable(t(tab13_7[ 1:10,]),col.names = rep("\\",10), caption="**Tabela 13.7**: Postos para o Exemplo 13.9"))
print(kable(t(tab13_7[11:20,]),col.names = rep("\\",10)))
print(kable(t(tab13_7[21:30,]),col.names = rep("\\",10)))
print(kable(t(tab13_7[31:40,]),col.names = rep("\\",10)))attach(tab13_7)
wilcox.test(media~tipo, alternative="less")A<-c(80,72,65,78,85)
B<-c(75,70,60,72,78)
operador<-1:5
tab13_8<-data.frame(operador,A,B)
t.test(A, B, alternative="greater", paired=TRUE, conf.level = 0.90)p.test(p1=168/400, p2=180/600,n1=400,n2=600,alpha=0.05, H1="p1!=p2")
IC_proporcoes(p1=168/400, p2=180/600,n1=400,n2=600,gama=0.95)
tab13_12A<-matrix(c(168,232,180,420),ncol=2,nrow=2,byrow = TRUE)
prop.test(x=tab13_12A, alternative = "two.sided")tab13_12<-data.frame(sujeito=1:26,rbind(
c(2.18 , 0.43),c(2.05 , 0.08),c(1.05 , 0.18),c(1.95 , 0.78),c(0.28 , 0.03),
c(2.63 , 0.23),c(1.5 , 0.2 ),c(0.45 , 0 ),c(0.7 , 0.05),c(1.3 , 0.3 ),
c(1.25 , 0.33),c(0.18 , 0 ),c(3.3 , 0.9 ),c(1.4 , 0.24),c(0.9 , 0.15),
c(0.58 , 0.1 ),c(2.5 , 0.33),c(2.25 , 0.33),c(1.53 , 0.53),c(1.43 , 0.43),
c(3.48 , 0.65),c(1.8 , 0.2 ),c(1.5 , 0.25),c(2.55 , 0.15),c(1.3 , 0.05),
c(2.65 , 0.25)))
names(tab13_12)[2:3]<-c("antes", "depois")
tab13_12$d<-tab13_12$antes-tab13_12$depois
tab13_12$postos<-rank(abs(tab13_12$d))
print(kable(tab13_12,caption="**Tabela 13.12**: Índices de placa bacteriana.",col.names = c("Sujeito","Antes(x_i)","Depois(y_i)","d_i=x_i-y_i","postos de |d_i|")))xp<-list(tab13_12$antes,tab13_12$depois)
boxplot(xp,pch="-", col="lightblue", border="black", boxwex=0.3, names=c("x","y"))diff_trats<-t.test(tab13_12$antes,tab13_12$depois, alternative="two.sided", paired=TRUE, conf.level = 0.95)
diff_tratsstripchart(tab13_12$d,method = "overplot", offset = 2, at=1,
pch = 1, col="darkred",ylab=NA,
xlim=c(0,3))
par(new=TRUE)
plot(x=c(0,mean(tab13_12$d)),
y=c(0.8,0.8),xlim=c(0,3),pch=c(20,20),ylim=c(0,2),ylab="",xlab="", col="darkblue")
par(new=TRUE)
plot(x=c(diff_trats$conf.int[1],diff_trats$conf.int[2]), col="blue",
y=c(0.8,0.8),xlim=c(0,3),pch=c("|","|"),ylim=c(0,2), type="b",ylab="",xlab="")
text(x=c(0,mean(tab13_12$d)),y=c(0.7,0.7), labels=c("H_0","Média(d)"), col="lightblue4")boxplot(tab13_12$d,pch="|", col="lightblue", border="black", boxwex=0.3,horizontal = TRUE )
par(new=TRUE)
plot(x=c(0,mean(tab13_12$d)),
y=c(0.8,0.8),xlim=c(0,3),pch=c(20,20),ylim=c(0,2),ylab="",xlab="", col="darkblue",xaxt='n', yaxt='n')
par(new=TRUE)
plot(x=c(diff_trats$conf.int[1],diff_trats$conf.int[2]), col="blue",
y=c(0.8,0.8),xlim=c(0,3),pch=c("|","|"),ylim=c(0,2), type="b",ylab="",xlab="",
xaxt='n', yaxt='n')
text(x=c(0,mean(tab13_12$d)),y=c(0.7,0.7), labels=c("H_0","Média(d)"), col="lightblue4")f_obs = c(43,49,56,45,66,41)
f_esp = c(50,50,50,50,50,50)
q_obs = sum((f_obs-f_esp)^2)/50
q_c=qchisq(0.95,5)
if(q_obs > q_c){
print("Rejeito H0")
}else print("Aceito H0")
chisq.test(f_obs)c_hum<-c(15,20,30,20,15)
c_bio<-c( 8,23,18,34,17)
tab14_2<-rbind(c_hum,c_bio)
test_tab14_2=chisq.test(tab14_2)
test_tab14_2
tab14_8 = rbind(test_tab14_2$expected,
total=apply(test_tab14_2$expected,2,sum))
tab14_8=rbind(test_tab14_2$expected, total=apply(test_tab14_2$expected,2,sum))
print(kable(tab14_8, caption = "**Tabela 14.8**: Frequências absolutas sob $H_0$", col.names = c("A","B","C","D","E")))M<-data.frame(uso_hospital=c("usaram_hospital", "nao_usaram_hospital"),homens=c(100,900),mulheres=c(150,850), row.names = 1)
Macidentes<-c(32,40,20,25,33)
test_tab14_5<-chisq.test(acidentes)
q_i<-(acidentes-test_tab14_5$expected)^2/test_tab14_5$expected
tab14_5<- rbind(acidentes, e_i=test_tab14_5$expected, q_i)
tab14_5<-cbind(tab14_5,c(sum(tab14_5[1,]),sum(tab14_5[2,]),sum(tab14_5[3,])))
test_tab14_5x<-0:10
n_k<- c(57,203,383,525,532,408,273,139,45,27,16)
np_k<-c(54.399,210.523,407.361,525.496,508.418,393.515,253.817,140.325,67.882,29.189,17.075)
tab6_13<-cbind(x,n_k,np_k)
teste_quiquad(n_k,np_k,0.05,10)dados.ex14_6<-c( 1.04, 1.73, 3.93, 4.44, 6.37, 6.51,
7.61, 7.64, 8.18, 8.48, 8.57, 8.65,
9.71, 9.87, 9.95, 10.01, 10.52, 10.69,
11.72, 12.17, 12.61, 12.98, 13.03, 13.16,
14.11, 14.60, 14.64, 14.75, 16.68, 22.14)
classe<-numeric()
classe[1]<-sum(dados.ex14_6<6.63)
classe[2]<-sum(dados.ex14_6 > 6.63 & dados.ex14_6 < 10)
classe[3]<-sum(dados.ex14_6 > 10 & dados.ex14_6 < 13.37)
classe[4]<-sum(dados.ex14_6 > 13.37)
test_tab14_6<- chisq.test(classe)tab14_6<-cbind(rbind(classe,e_i=test_tab14_6$expected),c(sum(classe),sum(test_tab14_6$expected)))
teste_quiquad(classe,test_tab14_6$expected,0.05,3)
chisq.test(classe)par(mfrow=c(1,2))
hist(dados.ex14_6,col="lightblue3",border = "white",xlab="Dados Normal")
boxplot(dados.ex14_6, xlab="Dados Normal",col="lightblue", pch=16, border="darkgray")tab14_10<-rbind(P_1T=c(29,60,9,2),P_2C=c(37,44,13,6))
chisq.test(tab14_10)# Dados das tabelas 4.8 e 4.9
dados.tab4_8<-data.frame(rbind(
matrix(rep(c("1.Consumidor","1.São Paulo"),times=214),ncol=2,byrow=T),
matrix(rep(c("1.Consumidor","2.Paraná"),times=51),ncol=2,byrow=T),
matrix(rep(c("1.Consumidor","3.Rio G. do Sul"),times=111),ncol=2,byrow=T),
matrix(rep(c("2.Produtor","1.São Paulo"),times=237),ncol=2,byrow=T),
matrix(rep(c("2.Produtor","2.Paraná"),times=102),ncol=2,byrow=T),
matrix(rep(c("2.Produtor","3.Rio G. do Sul"),times=304),ncol=2,byrow=T),
matrix(rep(c("3.Escola","1.São Paulo"),times=78),ncol=2,byrow=T),
matrix(rep(c("3.Escola","2.Paraná"),times=126),ncol=2,byrow=T),
matrix(rep(c("3.Escola","3.Rio G. do Sul"),times=139),ncol=2,byrow=T),
matrix(rep(c("4.Outras","1.São Paulo"),times=119),ncol=2,byrow=T),
matrix(rep(c("4.Outras","2.Paraná"),times=22),ncol=2,byrow=T),
matrix(rep(c("4.Outras","3.Rio G. do Sul"),times=48),ncol=2,byrow=T)))
colnames(dados.tab4_8)<-c("tipo_de_cooperativa","estado")
attach(dados.tab4_8)CrossTable(estado,tipo_de_cooperativa,
prop.r=FALSE, prop.c=FALSE, prop.t=FALSE, prop.chisq=FALSE, expected=TRUE,
digits=0)anos_exp<-c(2,4,5,6,8)
n_clientes<-c(48,56,64,60,72)
cor.test(anos_exp,n_clientes,alternative="two.sided", method="pearson",conf.level=0.95)ks.test(dados.ex14_6,y="pnorm",mean=10,sd=5)
# Gráfico qxq da distribuição emprírica e da distribuição Normal(10,25)
qqnorm(dados.ex14_6,cex=2,pch=20,col="darkblue", xlab="Quantis da Normal Padrão", ylab="Quantis dos dados", main="")
qqline(dados.ex14_6)attach(tab15_1)
summary2(tempo_Y[sexo_W=="H"])
summary2(tempo_Y[sexo_W=="M"])
summary_sexo<-sapply(levels(sexo_W),function(sex){summary2(tempo_Y[sexo_W==sex])})
summary_sexo
summary_idade<-sapply(levels(idade_X),function(age){summary2(tempo_Y[idade_X==age])})
summary_idademean_idade<-rep(0,20)
for(age in levels(idade_X)){
mean_idade[mean_idade==0]<-((idade_X==age)*mean(tempo_Y[idade_X==age]))[mean_idade==0]
}
mean_sexo<-rep(0,length(sexo_W))
for(sex in levels(sexo_W)){
mean_sexo[mean_sexo==0]<-((sexo_W==sex)*mean(tempo_Y[sexo_W==sex]))[mean_sexo==0]
}
e1<-tempo_Y-mean(tempo_Y)
e2<-tempo_Y-mean_sexo
e3<-tempo_Y-mean_idade
tab15_4<-cbind(tab15_1[,1:4],e1,e2,e3)IC_t( mean(tempo_Y[sexo_W=="H"]),mean(tempo_Y[sexo_W=="M"]),
var(tempo_Y[sexo_W=="H"]),var(tempo_Y[sexo_W=="M"]),10,10,gama=0.95)
IC_diff_medias(mean(tempo_Y[sexo_W=="M"]),mean(tempo_Y[sexo_W=="H"]),
var(tempo_Y[sexo_W=="M"]),var(tempo_Y[sexo_W=="H"]),10,10,gama=0.95)fit15_4_sexo<-anova(lm(tempo_Y~sexo_W))
par(mfrow=c(2,2))
plot(aov(lm(tempo_Y~sexo_W)))qf(0.95,1,18)
if(fit15_4_sexo$`F value`[1]>qf(0.95,1,18)){
print("Rejeito H_0")
} else print( "Aceito H0")boxplot(tempo_Y~idade_X, col="lightblue3")fit15_4_idade<-anova(lm(tempo_Y~idade_X))
par(mfrow=c(2,2))
plot(aov(lm(tempo_Y~idade_X)))plot(as.vector(idade_X),e3,col="darkred",pch=20,xaxt="n")
axis(1,at=c(19.6,24.6,29.6,34.6,39.6),labels=c("20 anos","25 anos","30 anos","35 anos","40 anos"))
abline(v=c(19.6,24.6,29.6,34.6,39.6), col="darkgrey", lty=2)qf(0.95,4,15)
if(fit15_4_idade$`F value`[1]>qf(0.95,4,15)){
print("Rejeito H_0")
} else print( "Aceito H0")medias_idade<-sapply(levels(idade_X),function(age){mean(tempo_Y[idade_X==age])})
medias_idade
diff(medias_idade)v<-levels(idade_X)
SQ<-fit15_4_idade$`Sum Sq`[1] #fit15_4_idade ;e o modelo ajustado e $`Sum Sq`[1] é a SQ para a idade.
ics_diff_idades<-matrix(NA,10,2)
cont=1
lab_ics<-character(0)
vec_diff_idades<-numeric()
for( i in 1:4)
for(j in (i+1):5){
ics_diff_idades[cont,]<-IC_t_diff(mean(tempo_Y[idade_X==v[i]]),mean(tempo_Y[idade_X==v[j]]),SQ,20,4,4,gama = 0.995)
vec_diff_idades[cont]<-mean(tempo_Y[idade_X==v[i]])-mean(tempo_Y[idade_X==v[j]])
lab_ics[cont]<-paste(v[i],v[j],sep="_")
cont=cont+1
}
rownames(ics_diff_idades)<-lab_ics
ics_diff_idades
### Grafico com os ICs
plot(1:10,vec_diff_idades, ylim=c(min(ics_diff_idades)-1,max(ics_diff_idades)+1), type="p",xaxt="n",xlab="pares",ylab="ICs")
axis(1,at=1:10,labels=lab_ics) # labels no eixo x
for (i in 1:10) {
if(ics_diff_idades[i,1]>0 || ics_diff_idades[i,2]<0) {color="red"} else color="darkblue"
lines(x =c(i,i), y=t(ics_diff_idades)[,i], col=color) # Plotando os ICs
}bartlett.test(x=tempo_Y, g=idade_X)fit15_4_idade
ics_idades<-matrix(NA,5,2)
cont=1
for(age in levels(idade_X)){
ics_idades[cont,]<-IC(n = 4, mean(tempo_Y[idade_X==age]),sd(tempo_Y[idade_X==age]),gama=0.95)
cont=cont+1
}
print(kable(cbind(medias_idade,ics_idades),caption="Média e Intervalos de confiança por grupo de idade",col.names = c("Média","Limite Inferior", "Limite Superior")))
### gráfico dos ICs
plot(c(20,25,30,35,40),medias_idade, ylim=c(min(ics_idades)-1,max(ics_idades)+1), type="p",xlab="Idade",ylab="ICs")
for (i in 1:5) lines(x =c(v[i],v[i]), y=t(ics_idades)[,i]) # Plotando os ICsboxplot(e3~as.vector(idade_X),col="lightblue3",pch=20)
boxplot(e3,col = "lightblue3")attach(tab15_1)
n_idade_X<-as.numeric(as.character(idade_X))
plot(n_idade_X,tempo_Y, pch=20, xlab="Idade(x)", ylab="Estímulo(y)", col="darkblue")
abline(lm(tempo_Y~n_idade_X), lwd=2, col="red")fit16_1<-lm(tempo_Y~n_idade_X)
fit16_1
# fit16_1$coefficients nos dá os coeficientes ajustados, sendo fit16_1$coefficients[1] o intercepto e fit16_1$coefficients [2] a inclinação.
y_20=fit16_1$coefficients[1]+fit16_1$coefficients[2]*20
y_20
predict(object = fit16_1)
y_33=fit16_1$coefficients[1]+fit16_1$coefficients[2]*33
y_33
# Na função `predict` é possivel também calcular os valores preditos para valores específicos das variáveis explicativas:
newdata=data.frame(n_idade_X=20) # `data.frame` contém os dados das variáveis explicativas utilizados para calcular as previsões.
predict(fit16_1, newdata)tab16_1<-cbind(tab15_1[,1:3],resid_fit16_1=resid(fit16_1))
print(kable(tab16_1,caption="**Tabela 16.1**: Resíduos para o modelo (16.18).", col.names = c("Tempo de Reação", "Sexo", "Idade", "Resíduos")))
anova_fit16_1<-anova(fit16_1)
SQRes=sum(tab16_1$resid_fit16_1^2)
S2e=SQRes/(20-2)
SQRes
S2e
sqrt(S2e)
2*sqrt(S2e)s2<-var(tempo_Y)
s2summary(fit16_1)$r.squared anova(fit16_1)confint(fit16_1)summary(fit16_1)newdata=data.frame(n_idade_X=28) # `data.frame` contém os dados das variáveis explicativas utilizados para calcular as previsões.
predict(fit16_1, newdata, interval="predict") attach(tab16_1)
nui=1/length(n_idade_X)+(n_idade_X-mean(n_idade_X))^2/sum((n_idade_X-mean(n_idade_X))^2)
ri=resid_fit16_1/(sqrt(S2e)*sqrt(1-nui))
tab16_4<-cbind(tab16_1[3:4],zi=resid_fit16_1/sqrt(S2e),ri)par(mfrow=c(1,2))
plot(n_idade_X,resid_fit16_1, cex=2, col="darkblue", pch=20, xlab="Idade", ylab="Resíduos",main="(a)")
abline(h=0)
plot(n_idade_X,tab16_4$zi, cex=2, col="darkblue", pch=20, xlab="Idade", ylab="Resíduos padronizados",main="(b)")
abline(h=0)hist(resid_fit16_1, xlab="Resíduos", ylab="", col="lightblue3", border="white", main="")qqnorm(resid_fit16_1, cex=1.5, col="darkblue", xlab="Quantis da normal padrão", ylab="Quantis dos resíduos", pch=20)
qqline(resid_fit16_1, col="darkred")Análise de Resíduos no R:
par(mfrow=c(2,2))
plot(fit16_1)fit16_5<-lm(tab16_5$yi~tab16_5$xi)
fit16_5attach(tab16_5)
plot(xi,yi, pch=20, xlab="Temperatura", ylab="Vapor", col="darkblue",main="(a)")
abline(lm(yi~xi), lwd=2, col="red")par(mfrow=c(1,2))
plot(xi,resid(fit16_5), cex=2, col="darkblue", pch=20, xlab="Temperatura", ylab="Resíduos",main="(b)")
abline(h=0)qqnorm(resid(fit16_5), cex=1.5, col="darkblue", xlab="Quantis da normal padrão", ylab="Quantis dos resíduos", pch=20,main="(c)")
qqline(resid(fit16_5), col="darkred")#Ajuste do modelo para os dados do Exemplo 16.9 sem o intercepto é ajustado adicionando `-1` à equação:
fit16_9<-lm(tab16_9A$y~tab16_9A$x-1)
fit16_9
summary(fit16_9)
S2e=sum(resid(fit16_9)^2)/(9-1)
S2e
sqrt(S2e)
confint(fit16_9)
plot(tab16_9A$x,tab16_9A$y, pch=20, xlab="x", ylab="y", col="darkblue")
abline(lm(tab16_9A$y~tab16_9A$x), lwd=2, col="red")ano=seq(1961,1979,2)
inflacao_Y=c(9,24,72,128,192,277,373,613,1236,2639)
tab16_6<-data.frame(
ano=ano,
t=ano-mean(ano),
inflacao_Y=inflacao_Y,
log_Y=log(inflacao_Y))plot(t,inflacao_Y, pch=20, col="darkblue", ylim=c(0,max(inflacao_Y)+5),ylab="Inflação",xlab="Ano(t)")
par(new=T)
plot(t,exp(predict(fit16_10)),pch="+", col="darkred", ylim=c(0,max(inflacao_Y)+5),ylab="",xlab="")fit16_10<-lm(log_Y~t)
fit16_10
exp(fit16_10$coefficients[1])plot(tab16_6$t,tab16_6$log_Y, pch=20, xlab="x", ylab="y", col="darkblue")
abline(lm(tab16_6$log_Y~tab16_6$t), lwd=2, col="red")
resid_exp=inflacao_Y-exp(predict(fit16_10))
tab16_7<-data.frame(t=tab16_6$t,resid_reta=resid(fit16_10),resid_exp)
print(kable(tab16_7,caption="**Tabela 16.7**: Resíduos para os modelos linear e exponencial.", col.names=c("t", "Resíduos Reta", "Resíduos Exponencial")))plot(tab16_7$t,tab16_7$resid_reta, xlab="ano", ylab="Resíduos reta", pch=20)
abline(h=0,lty=2, col="gray")plot(tab16_7$t,tab16_7$resid_exp, xlab="ano", ylab="Resíduos exponencial", pch=20)
abline(h=0,lty=2, col="gray")par(mfrow=c(1,2))
hist(tab16_7$resid_reta,col="lightblue3",border="white", main="(a)",xlab="",ylab="Resíduos reta")
hist(tab16_7$resid_exp,col="lightblue3",border="white", main="(b)",xlab="",ylab="Resíduos exponencial")par(mfrow=c(1,2))
qqnorm(tab16_7$resid_reta,col="darkblue", main="(a)",xlab="Quantis da normal padrão",ylab="Resíduos reta", pch=20,cex=1.5)
qqline(tab16_7$resid_reta, col="darkred")
qqnorm(tab16_7$resid_exp,col="darkblue", main="(b)",xlab="Quantis da normal padrão",ylab="Resíduos exponencial", pch=20,cex=1.5)
qqline(tab16_7$resid_exp, col="darkred")# Definin
grupo_E=c(2,1,4,3,5,8,6)
grupo_C=c(7,12,10,11,9,14)
grupo_D=c(16,13,15,18,17,20,19)
fit16_11<-reg_resist(data.frame(x=n_idade_X,y=tab15_1$tempo_Y))
fit16_11cd_mercado <- read.table("cd-mercado.csv",h=T,skip=4,sep=";",dec=",") # Leitura dos dados
attach(cd_mercado)
fit16_12<-lm(telebras[1:39]~indice[1:39])
plot(indice[1:39],telebras[1:39],ylab="Tel",xlab="IBV",pch=16,col="darkblue")
abline(fit16_12)summary(fit16_12)
anova(fit16_12)resid_fit16_12<-resid(fit16_12)
par(mfrow=c(2,2))
qqnorm(resid_fit16_12, pch=20, col="darkblue")
qqline(resid_fit16_12, col="darkred")
plot.ts(resid_fit16_12, pch=20, col="darkblue")
abline(h=c(-0.7932,0.7932), col="red", lty=2)
hist(resid_fit16_12, col="lightblue3", border = "white")
plot(fitted(fit16_12),resid_fit16_12, pch=20, col="darkblue")
abline(h=0)
par(mfrow=c(2,2))
plot(fit16_12)attach(tab16_8)
fit_16_13<-lm(vt~t)
fit_16_13
fit_resist16_13<-reg_resist(data.frame(x=t,y=vt))
fit_resist16_13