Métodos Computacionais para Estatística 2
Github com os arquivos
)
)
Lista 01
library(readxl)
library(magrittr)
library(readr)
library(haven)
library(stringr)
library(dplyr)Questão 1
(a)
(8/22)^(-3) + (gamma(4)/22)## [1] 21.0696(b)
sqrt(25/102) + log10(33)## [1] 2.013588(c)
(100/22)^(-1/4) + exp(-2/3) - 2/3## [1] 0.5316166(d)
choose(10, 4) - factorial(7)## [1] -4830(e)
beta(2, 3) + gamma(9) + abs(-5) - log(22)## [1] 40321.99Questão 2
(a)
seq(1, 10)## [1] 1 2 3 4 5 6 7 8 9 10(b)
seq(1, 100, 3)## [1] 1 4 7 10 13 16 19 22 25 28 31 34 37 40 43 46 49
## [18] 52 55 58 61 64 67 70 73 76 79 82 85 88 91 94 97 100(c)
seq(0, 1000, 50)## [1] 0 50 100 150 200 250 300 350 400 450 500 550 600 650
## [15] 700 750 800 850 900 950 1000(d)
rep(c(1, 22, 5, 3), c(6, 2, 2, 4))## [1] 1 1 1 1 1 1 22 22 5 5 3 3 3 3(e)
rep(c(9, 10, 11, 12), each=2, times=3)## [1] 9 9 10 10 11 11 12 12 9 9 10 10 11 11 12 12 9 9 10 10 11 11 12
## [24] 12Questão 3
A <- matrix(data = c(
rep(c(0,1),500),
1:1000,
rep(c(10,20,30,40,50), 200)
), ncol = 3)(a)
A[seq(1,1000,2), c(1,3)] %>% head()## [,1] [,2]
## [1,] 0 10
## [2,] 0 30
## [3,] 0 50
## [4,] 0 20
## [5,] 0 40
## [6,] 0 10(b)
A[,1] <- rep(1:200, each=5)(c)
A[seq(0,1000,2),] <- 0Questão 4
base_Q4_1 <- tibble(
Nome = c('Joao','Maria','Ana','Fabia','Rodrigo','Renato','Luciana','Guilherme','Gabriel','Diogo'),
Idade = c(27,22,21,37,29,27,21,18,19,25),
Altura = c(1.83,1.53,1.79,1.58,1.65,1.70,1.51,1.66,1.72,1.83),
Genero = c("M","F","F","F","M","M","F","M","M","M"),
Estado_Civil = c("Solteiro","Solteiro","Casado","Solteiro","Casado","Solteiro","Solteiro","Solteiro","Casado","Casado"),
Bairro = c('Icarai','Inga','Botafogo','Lagoa','Boa Viagem','Leblon','Leblon','Inga','Icarai','Botafogo')
)
base_Q4_2 <- tibble(
nome = c('Joao','Maria','Ana','Fabia','Rodrigo','Renato','Gabriel','Diogo'),
peso = c(40,65,77,63,78,80,83,77)
)
base_Q4_3 <- tibble(
Nome = c('Joao','Maria','Ana','Fabia','Rodrigo','Renato','Gabriel','Diogo','Vicente','Fernando'),
opniao = c("Contra","Contra","A favor","Contra","A favor","Contra","A favor","A favor","A favor","A favor")
)
base_Q4_4 <- tibble(
bairro = c('Inga','Icarai','Boa Viagem','Botafogo','Leblon','Copacabana','Ipanema','Lagoa','Gavea','Sao Francisco'),
Cidade = c("Niterói","Niterói","Niterói","Rio de Janeiro","Rio de Janeiro","Rio de Janeiro","Rio de Janeiro","Rio de Janeiro","Rio de Janeiro","Rio de Janeiro")
)(a)
base_Q4_a <- inner_join(base_Q4_1, base_Q4_2, by=c("Nome"="nome"))(b)
saveRDS(base_Q4_a, "Lista 01/base_Q4_a.rds")(c)
base_Q4_c <- right_join(base_Q4_1, base_Q4_2, by=c("Nome"="nome"))(d)
write_csv(base_Q4_c, "Lista 01/base_Q4_c.csv")(e)
base_Q4_e <- full_join(base_Q4_1, base_Q4_2, by=c("Nome"="nome")) %>%
full_join(base_Q4_3, by="Nome")(f)
write_dta(base_Q4_e, "Lista 01/base_Q4_e.dta")(g)
base_Q4_g <- left_join(base_Q4_1, base_Q4_4, by=c("Bairro"="bairro"))(h)
write_delim(base_Q4_g, "Lista 01/base_Q4_g.txt")(i)
# file.size(dir("Lista 01")[1]) %>% min()Questão 5
Banco1 <- read_excel("Lista 01/Banco1.xls")(a)
Banco1$Altura %>% max()## [1] 1.98(b)
Banco1$Peso[which.max(Banco1$Altura)]## [1] 95(c)
Banco1$Grupo[which.min(Banco1$Peso)]## [1] 1(d)
Banco1$Grupo %<>% factor(labels = c("Placebo", "Tratamento A"))(e)
write_csv(Banco1, "Lista 01/Banco1_mod.csv")Questão 6
dados_pop <- read_csv("Lista 01/populacaototaljovem2010.csv")
dados2010 <- read_excel("Lista 01/dados2010.xls")(a)
dados2010 %>% distinct(codmun)## # A tibble: 1,880 x 1
## codmun
## <dbl>
## 1 110000
## 2 110001
## 3 110002
## 4 110003
## 5 110004
## 6 110005
## 7 110006
## 8 110007
## 9 110008
## 10 110009
## # ... with 1,870 more rows(b)
dados2010 %>%
group_by(codmun) %>%
summarise(tot_estupros = sum(Estupros, na.rm = T))## # A tibble: 1,880 x 2
## codmun tot_estupros
## <dbl> <dbl>
## 1 110000 0
## 2 110001 8
## 3 110002 50
## 4 110003 0
## 5 110004 25
## 6 110005 7
## 7 110006 0
## 8 110007 1
## 9 110008 6
## 10 110009 2
## # ... with 1,870 more rows(c)
dados2010 %<>%
mutate(idade.cat = cut(Idade %>% as.integer(), breaks = c(-Inf, 18, Inf), labels = c(1, 2), right=F))
dados2010 %>%
group_by(codmun, idade.cat) %>%
summarise(tot_afogamentos = sum(Afogamentos, na.rm = T))## # A tibble: 3,365 x 3
## # Groups: codmun [?]
## codmun idade.cat tot_afogamentos
## <dbl> <fct> <dbl>
## 1 110000 1 1
## 2 110001 1 1
## 3 110001 2 6
## 4 110002 1 5
## 5 110002 2 28
## 6 110003 1 0
## 7 110003 2 0
## 8 110003 <NA> 0
## 9 110004 1 3
## 10 110004 2 23
## # ... with 3,355 more rows(d)
dados2010 %>%
group_by(codmun) %>%
summarise(pop_total = length(Idade),
pop_jovem = sum(idade.cat==1, na.rm = T))## # A tibble: 1,880 x 3
## codmun pop_total pop_jovem
## <dbl> <int> <int>
## 1 110000 2 2
## 2 110001 47 3
## 3 110002 88 10
## 4 110003 10 1
## 5 110004 87 8
## 6 110005 40 3
## 7 110006 36 4
## 8 110007 13 2
## 9 110008 26 1
## 10 110009 49 2
## # ... with 1,870 more rows(e)
dados2010 %<>% rename("maior.idade" = "idade.cat")(f)
dados2010 %>%
group_by(codmun) %>%
filter(Idade<18)## # A tibble: 4,518 x 6
## # Groups: codmun [1,387]
## codmun Idade Estupros Espancamentos Afogamentos maior.idade
## <dbl> <chr> <dbl> <dbl> <dbl> <fct>
## 1 110000 14 0 0 1 1
## 2 110000 17 0 0 0 1
## 3 110001 0 0 0 0 1
## 4 110001 1 0 0 1 1
## 5 110001 14 0 0 0 1
## 6 110002 0 0 0 1 1
## 7 110002 1 0 0 1 1
## 8 110002 10 0 0 0 1
## 9 110002 11 0 0 0 1
## 10 110002 15 0 0 1 1
## # ... with 4,508 more rowsQuestão 7
base_mae <- read_dta("Lista 01/basemae.dta")(a)
base_mae %<>%
mutate(UF = str_sub(munic_res, end = 2),
GR = str_sub(munic_res, end = 1))(b)
base_mae %>%
group_by(UF) %>%
filter(UF == 21, ano >= 2000, ano <= 2005) %>%
summarise(taxa_media_estupro = mean(txestupro, na.rm = T),
taxa_media_agressao = mean(txagressao, na.rm = T))## # A tibble: 1 x 3
## UF taxa_media_estupro taxa_media_agressao
## <chr> <dbl> <dbl>
## 1 21 0.00313 5.45(c)
base_mae %>%
group_by(meso_res, meso_nome_res) %>%
summarise(num_mun = sum(length(munic_res), na.rm = T))## # A tibble: 138 x 3
## # Groups: meso_res [?]
## meso_res meso_nome_res num_mun
## <chr> <chr> <int>
## 1 "" "" 108
## 2 1101 Madeira-Guaporé 110
## 3 1102 Leste Rondoniense 462
## 4 1201 Vale do Juruá 88
## 5 1202 Vale do Acre 154
## 6 1301 Norte Amazonense 66
## 7 1302 Sudoeste Amazonense 176
## 8 1303 Centro Amazonense 330
## 9 1304 Sul Amazonense 110
## 10 1401 Norte de Roraima 88
## # ... with 128 more rows(d)
base_mae %>%
group_by(ano) %>%
summarise(taxa_media_suicidio = mean(txsuicidio, na.rm = T),
taxa_media_homicidio = mean(txhomicidio, na.rm = T))## # A tibble: 11 x 3
## ano taxa_media_suicidio taxa_media_homicidio
## <dbl> <dbl> <dbl>
## 1 2000 5.45 12.1
## 2 2001 6.04 13.1
## 3 2002 6.12 13.8
## 4 2003 6.00 14.3
## 5 2004 6.21 14.3
## 6 2005 6.40 14.6
## 7 2006 6.21 14.6
## 8 2007 6.48 15.1
## 9 2008 6.41 15.5
## 10 2009 6.36 16.1
## 11 2010 6.38 16.3(e)
base_mae %>%
group_by(ano, GR) %>%
filter(!is.na(txsuicidio), !is.na(txhomicidio)) %>%
summarise(taxa_media_suicidio = mean(txsuicidio, na.rm = T),
taxa_media_homicidio = mean(txhomicidio, na.rm = T))## # A tibble: 55 x 4
## # Groups: ano [?]
## ano GR taxa_media_suicidio taxa_media_homicidio
## <dbl> <chr> <dbl> <dbl>
## 1 2000 1 2.48 14.0
## 2 2000 2 2.42 11.2
## 3 2000 3 4.51 11.9
## 4 2000 4 10.9 9.55
## 5 2000 5 7.62 21.5
## 6 2001 1 3.41 16.6
## 7 2001 2 3.12 12.0
## 8 2001 3 5.11 12.6
## 9 2001 4 11.3 10.3
## 10 2001 5 7.80 22.1
## # ... with 45 more rows(f)
base_mae %>%
group_by(micro_res, micro_nome_res) %>%
summarise(pop_media = mean(populacao, na.rm = T),
tot_homicidios = sum(homicidio, na.rm = T),
tot_homicidios_paf = sum(homicidio_paf, na.rm = T))## # A tibble: 559 x 5
## # Groups: micro_res [?]
## micro_res micro_nome_res pop_media tot_homicidios tot_homicidios_paf
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 "" "" 2654. 87 60
## 2 11001 Porto Velho 67530. 2876 1903
## 3 11002 Guajará-Mirim 22476. 170 96
## 4 11003 Ariquemes 22750. 989 674
## 5 11004 Ji-Paraná 28230. 793 513
## 6 11005 Alvorada D'Oeste 19198. 144 97
## 7 11006 Cacoal 25658. 490 305
## 8 11007 Vilhena 19512. 357 175
## 9 11008 Colorado do Oeste 11282 82 37
## 10 12001 Cruzeiro do Sul 23310. 148 26
## # ... with 549 more rowsLista 02
library(readxl)
library(magrittr)
library(janitor)
library(dplyr)
library(readr)
library(stringr)
library(ggplot2)
library(gridExtra)Questão 1
data1 <- read_excel("Lista 02/Tratamento.xls", na = "999") %>% clean_names()(a)
data1$grupo %<>% factor(labels=c("Adesão", "Não Adesão"))
data1$situacao_atual_no_trabalho %<>% factor(labels=c("Empregado", "Desempregado"))
data1$frequenta_algum_grupo_religioso %<>% factor(labels=c("Não","Sim"))
data1$escolaridade %<>% factor(labels=c("analfabeto", "ensino fundamental", "ensino medio", "ensino superior"))(b)
data1 %<>%
mutate(faixa_etaria = cut(data1$idade, breaks=c(0,18,60,150), labels=c("Jovem", "Adulto", "3ª Idade")))(c)
class(data1)## [1] "tbl_df" "tbl" "data.frame"(d)
summary(data1)## grupo situacao_atual_no_trabalho idade
## Adesão :30 Empregado :17 Min. :67.00
## Não Adesão:26 Desempregado:39 1st Qu.:73.00
## Median :76.50
## Mean :76.32
## 3rd Qu.:79.25
## Max. :91.00
##
## escolaridade ansiedade depressao
## analfabeto : 4 Min. :15.0 Min. :10.0
## ensino fundamental:34 1st Qu.:19.0 1st Qu.:16.0
## ensino medio :10 Median :21.0 Median :18.0
## ensino superior : 8 Mean :20.8 Mean :17.6
## 3rd Qu.:23.0 3rd Qu.:20.0
## Max. :26.0 Max. :23.0
## NA's :1 NA's :1
## frequenta_algum_grupo_religioso faixa_etaria
## Não :29 Jovem : 0
## Sim :24 Adulto : 0
## NA's: 3 3ª Idade:56
##
##
##
## (e)
write_delim(data1, "Lista 02/Tratamento_Q1.txt")
write_csv(data1, "Lista 02/Tratamento_Q1.csv")
write_rds(data1, "Lista 02/Tratamento_Q1.rds")
files <- c("Lista 02/Tratamento_Q1.txt",
"Lista 02/Tratamento_Q1.csv",
"Lista 02/Tratamento_Q1.rds")
for (i in 1:2) {
ifelse (((file.size(files[i])) < (file.size(files[i+1]))), j <- 1, j <- i+1)
}; cat("O arquivo com menor tamanho é ",
str_sub((files[j]), start = 10),
", pois ao salvar o arquivo com a extensão .rds o R salva-o com a compressão gzip, o que diminui o tamanho e o tempo necessário para a leitura do mesmo.")## O arquivo com menor tamanho é Tratamento_Q1.rds , pois ao salvar o arquivo com a extensão .rds o R salva-o com a compressão gzip, o que diminui o tamanho e o tempo necessário para a leitura do mesmo.Questão 2
data2 <- read_rds("Lista 02/Tratamento_Q1.rds")(a)
data2 %>% select(idade, ansiedade, depressao) %>% write_delim("Lista 02/Tratamento_Q2a.txt")(b)
data2 %>%
group_by(escolaridade) %>%
filter(escolaridade=="ensino fundamental", idade>50, frequenta_algum_grupo_religioso=="Não") %>%
write_csv("Lista 02/Tratamento_Q2b.csv")(c)
data2 %>% arrange(idade)## # A tibble: 56 x 8
## grupo situacao_atual_~ idade escolaridade ansiedade depressao
## <fct> <fct> <dbl> <fct> <dbl> <dbl>
## 1 Não ~ Empregado 67 ensino medio 23 18
## 2 Não ~ Desempregado 67 ensino fund~ 22 23
## 3 Não ~ Desempregado 67 ensino medio 25 17
## 4 Não ~ Desempregado 67 ensino fund~ 25 23
## 5 Não ~ Desempregado 67 ensino fund~ 15 19
## 6 Não ~ Empregado 68 ensino supe~ 25 22
## 7 Ades~ Desempregado 69 ensino fund~ 20 14
## 8 Ades~ Desempregado 69 ensino medio 22 17
## 9 Ades~ Desempregado 70 ensino fund~ 19 12
## 10 Ades~ Desempregado 71 ensino medio 24 20
## # ... with 46 more rows, and 2 more variables:
## # frequenta_algum_grupo_religioso <fct>, faixa_etaria <fct>(d)
data2 %>% rename("pont.ansiedade" = "ansiedade")## # A tibble: 56 x 8
## grupo situacao_atual_~ idade escolaridade pont.ansiedade depressao
## <fct> <fct> <dbl> <fct> <dbl> <dbl>
## 1 Ades~ Empregado 73 ensino fund~ 23 15
## 2 Ades~ Empregado 72 ensino fund~ 25 18
## 3 Ades~ Empregado 80 analfabeto 21 20
## 4 Ades~ Empregado 78 ensino fund~ 20 15
## 5 Ades~ Empregado 81 ensino supe~ 21 20
## 6 Ades~ Empregado 77 ensino fund~ 19 16
## 7 Ades~ Empregado 90 ensino medio 15 18
## 8 Ades~ Empregado 81 ensino fund~ 21 19
## 9 Ades~ Desempregado 70 ensino fund~ 19 12
## 10 Ades~ Desempregado 74 ensino fund~ 16 18
## # ... with 46 more rows, and 2 more variables:
## # frequenta_algum_grupo_religioso <fct>, faixa_etaria <fct>(e)
data2 %<>% mutate(pont.total = ansiedade+depressao)(f)
data2 %>%
group_by(grupo, situacao_atual_no_trabalho) %>%
summarise(med_ansiedade = mean(ansiedade, na.rm = T),
med_depressao = mean(depressao, na.rm = T))## # A tibble: 4 x 4
## # Groups: grupo [?]
## grupo situacao_atual_no_trabalho med_ansiedade med_depressao
## <fct> <fct> <dbl> <dbl>
## 1 Adesão Empregado 20.6 17.6
## 2 Adesão Desempregado 20.3 15.9
## 3 Não Adesão Empregado 20.7 19.6
## 4 Não Adesão Desempregado 21.7 18.6(g)
data2 %>%
group_by(grupo, escolaridade, frequenta_algum_grupo_religioso) %>%
summarise(maximo = max(pont.total, na.rm = T),
mediana = median(pont.total, na.rm = T),
"5o_decil" = quantile(pont.total, 0.5, na.rm = T))## # A tibble: 15 x 6
## # Groups: grupo, escolaridade [?]
## grupo escolaridade frequenta_algum_grup~ maximo mediana `5o_decil`
## <fct> <fct> <fct> <dbl> <dbl> <dbl>
## 1 Adesão analfabeto Não 28 28 28
## 2 Adesão analfabeto Sim 41 38 38
## 3 Adesão ensino fundame~ Não 42 32.5 32.5
## 4 Adesão ensino fundame~ Sim 43 35 35
## 5 Adesão ensino medio Não 44 41.5 41.5
## 6 Adesão ensino medio Sim 37 35 35
## 7 Adesão ensino superior Sim 47 42 42
## 8 Não Ad~ analfabeto Não 38 38 38
## 9 Não Ad~ ensino fundame~ Não 48 38 38
## 10 Não Ad~ ensino fundame~ Sim 45 40 40
## 11 Não Ad~ ensino fundame~ <NA> 45 45 45
## 12 Não Ad~ ensino medio Não 43 41 41
## 13 Não Ad~ ensino superior Não 37 36.5 36.5
## 14 Não Ad~ ensino superior Sim 47 47 47
## 15 Não Ad~ ensino superior <NA> 43 43 43Questão 3
(a)
grafa <-
data1 %>%
ggplot() +
geom_bar(aes(x=grupo), fill="#8cb8ff", color="black") +
theme_bw() +
labs(x="Grupo", y="", title="Grupos")(b)
grafb <-
data1 %>%
ggplot() +
geom_bar(aes(x=situacao_atual_no_trabalho), fill="#a570d1", color="black") +
theme_bw() +
labs(x="Situação", y="", title="Situação atual do trabalho")(c)
grafc <-
data1 %>%
ggplot() +
geom_boxplot(aes(x=grupo, y=idade, color=grupo)) +
theme_bw() +
labs(x="Grupo", y="Idade", title="Idade por Grupos") +
theme(legend.position = "none") ##ou guides(color=F)(d)
data1 %>%
filter(frequenta_algum_grupo_religioso=="Sim") %>%
ggplot() +
geom_boxplot(aes(x=ansiedade, y=depressao, group=frequenta_algum_grupo_religioso)) +
theme_bw() +
labs(x="Ansiedade", y="Depressão", title="Ansiedade e Depressão segundo frequenta algum grupo religioso")
(e)
data1 %>%
ggplot() +
geom_histogram(aes(x=idade), fill="#8cb8ff") +
theme_bw() +
labs(x="Idade", y="", title="Idade")
(f)
summary(data1)## grupo situacao_atual_no_trabalho idade
## Adesão :30 Empregado :17 Min. :67.00
## Não Adesão:26 Desempregado:39 1st Qu.:73.00
## Median :76.50
## Mean :76.32
## 3rd Qu.:79.25
## Max. :91.00
##
## escolaridade ansiedade depressao
## analfabeto : 4 Min. :15.0 Min. :10.0
## ensino fundamental:34 1st Qu.:19.0 1st Qu.:16.0
## ensino medio :10 Median :21.0 Median :18.0
## ensino superior : 8 Mean :20.8 Mean :17.6
## 3rd Qu.:23.0 3rd Qu.:20.0
## Max. :26.0 Max. :23.0
## NA's :1 NA's :1
## frequenta_algum_grupo_religioso faixa_etaria
## Não :29 Jovem : 0
## Sim :24 Adulto : 0
## NA's: 3 3ª Idade:56
##
##
##
## graf1 <-
data1 %>%
ggplot(aes(x=idade)) +
geom_histogram(aes(y=..density..), bins=10) +
stat_function(fun=dnorm, args=list(mean=mean(data1$idade), sd=sd(data1$idade)), color="red") +
theme_bw() +
labs(x="Idade", y="", title="Idade")
# Idade é uma variável quantitativa discreta
## Escolaridade é uma variável qualitativa, logo não possui densidade
graf2 <-
data1 %>%
ggplot(aes(x=depressao, na.rm=T)) +
geom_histogram(aes(y=..density..), bins=10) +
stat_function(fun=dnorm, args=list(mean=mean(na.omit(data1$depressao)), sd=sd(na.omit(data1$depressao))), color="red") +
theme_bw() +
labs(x="Depressão", y="", title="Depressão")
# Depressão é uma variável quantitativa discreta
graf3 <-
data1 %>%
na.omit(.$ansiedade) %>%
ggplot(aes(x=ansiedade, na.rm=T)) +
geom_histogram(aes(y=..density..), bins=10) +
stat_function(fun=dnorm, args=list(mean=mean(na.omit(data1$ansiedade)), sd=sd(na.omit(data1$ansiedade))), color="red") +
theme_bw() +
labs(x="Ansiedade", y="", title="Ansiedade")
# Ansiedade é uma variável quantitativa discreta
grid.arrange(graf1, graf2, graf3, ncol=2)
(g)
table(data1$situacao_atual_no_trabalho, data1$grupo) ##
## Adesão Não Adesão
## Empregado 8 9
## Desempregado 22 17table(data1$situacao_atual_no_trabalho, data1$grupo) %>% prop.table(2)##
## Adesão Não Adesão
## Empregado 0.2666667 0.3461538
## Desempregado 0.7333333 0.6538462data1 %>%
ggplot(aes(situacao_atual_no_trabalho)) +
geom_bar(aes(fill=grupo), position = "fill") +
theme_minimal() +
labs(x="Situação atual do trabalho", y="")
(h)
grid.arrange(grafa, grafb, grafc, ncol=2)
Questão 4
data1$depressao %>% summary()## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 10.0 16.0 18.0 17.6 20.0 23.0 1a4 <- data1$depressao %>% na.omit()
data1 %>%
ggplot() +
geom_boxplot(aes(y=depressao)) +
scale_y_continuous(breaks = seq(min(a4), max(a4), 2)) +
theme_bw()
data1$depressao %>% boxplot()
data1$frequenta_algum_grupo_religioso %>% summary()## Não Sim NA's
## 29 24 3data1 %>% filter(!is.na(frequenta_algum_grupo_religioso)) %>%
ggplot() +
geom_bar(aes(x=frequenta_algum_grupo_religioso, fill=frequenta_algum_grupo_religioso)) +
labs(x="Frequenta algum grupo religioso", y="", title="Frequenta algum grupo religioso") +
theme_bw() +
guides(fill=F)
Questão 5
(a)
data1 %<>% mutate(cat_idade = cut(idade, breaks=c(-Inf,70,90,+Inf), labels=c(1,2,3)))(b)
data1 %<>% mutate(cat_depre = cut(depressao, breaks=c(0,10,20,30), labels=c(1,2,3)))Questão 6
(a)
a1 <- function(x, lambda) {
lambda*exp(-lambda*x) * I(x>=0)
}(b)
a2 <- function(x) {
x^2 * I(x>-5 & x<5)
}(c)
a3 <- function(x) {
(1/10)*exp(-1)*(x-10) * I(x>2 & x<18)
}(d)
a4 <- function(alfa, beta, x) {
((beta^alfa)/gamma(alfa)) * x^(alfa-1) * exp(-beta * x) * I(x>=0)
}(e)
a5 <- function(alfa, beta, x) {
(gamma(alfa+beta)/((gamma(alfa)*gamma(beta)))) * x^(alfa-1) * (1-x)^(beta-1) * I(x<1 & x>0)
}
# a5 <- function(a, b, x) {
# ((x^(a-1))*((1-x)^(b-1)))/beta(a,b) * I(x<1 & x>0)
# }Questão 7
ggplot(data.frame(x=c(0,10)), aes(x=x)) +
stat_function(fun = a1, args = list(lambda=1)) +
labs(y="", title="Destribuição exponencial com lambda=1") +
theme_minimal()
ggplot(data.frame(x=c(-5,5)), aes(x=x)) +
stat_function(fun = a2) +
labs(y="f(x)", title="Função x^2 no intervalo (-5,5)") +
theme_minimal()
ggplot(data.frame(x=c(2,18)), aes(x=x)) +
stat_function(fun = a3) +
labs(y="f(x)") +
theme_minimal()
ggplot(data.frame(x=c(0,30)), aes(x=x)) +
stat_function(fun = a4, args = list(alfa=2, beta=1)) +
labs(y="", title="Distribuição Gama(2,1)") +
theme_minimal()
ggplot(data.frame(x=c(0,1)), aes(x=x)) +
stat_function(fun = a5, args = list(a=7, b=2)) +
labs(y="", title="Distribuição Beta(7,2)") +
theme_minimal()
Questão 8
dev_pad <- function(x) {
sqrt(var(x))
}
vet <- c(1:5)
sd(vet)## [1] 1.581139dev_pad(vet)## [1] 1.581139Questão 9
f9 <- function(x) {
cat(paste0("Soma = ", sum(x), ". Média = ", mean(x), "."))
}
f9(vet)## Soma = 15. Média = 3.Questão 10
IMC <- function(peso, altura) {
peso/((altura)^2)
}
IMC(80,1.80)## [1] 24.69136Questão 11
IMC <- function(peso, altura) {
imc <- peso/((altura)^2)
if (imc>=0 && imc<=18.5) {
cat("IMC: ", imc %>% round(), "- Desnutrido")
}else{
if(imc>18.5 && imc<=25){
cat("IMC: ", imc %>% round(), "- Normal")
}else{
cat("IMC: ", imc %>% round(), "- Obeso")
}
}
}
IMC(80,1.80)## IMC: 25 - NormalQuestão 12
(a)
{
soma=1
for (i in seq(32,20*100,20)) {
soma=soma+(1/i)
}
soma
}## [1] 1.224006(b)
{
soma=1
vet <- seq(11,10*100,10)
for (j in seq(1,50,2)) {
soma=soma-(1/vet[j])+(1/vet[j+1])
}
soma
}## [1] 0.9390823(c)
{
soma=0
for (i in 1:100) {
soma=soma+(i/15)
}
soma
}## [1] 336.6667(d)
{
soma=1/10
for (i in seq(10,10*99,10)) {
soma=soma+(i/10)
}
soma
}## [1] 4950.1(e)
{
soma=0
for (i in 1:100) {
soma=soma+sqrt(log(i))
}
soma
}## [1] 188.104Questão 13
(a)
(function(x) x^3) %>%
integrate(lower = 0, upper = 10)## 2500 with absolute error < 2.8e-11(b)
fb <- function(x) (3/5)*(x^3 + x)
integrate(fb, lower = 0, upper = 5)[[1]] + integrate(fb, lower = 7, upper = 10)[[1]] + integrate(fb, lower = 11, upper = 15)[[1]]## [1] 6685.2(c)
(function(x) x^12 * (1-8)^8) %>%
integrate(lower = 0, upper = 1)## 443446.2 with absolute error < 4.9e-09(d)
(function(x) 3*exp(-3*x)) %>%
integrate(lower = 0, upper = 100)## 1 with absolute error < 4.2e-08Questão 14
(a)
quote(x^3) %>% D("x")## 3 * x^2(b)
quote(cos(2*x)+exp(-3*x)) %>% D("x")## -(exp(-3 * x) * 3 + sin(2 * x) * 2)(c)
quote(3*x+log(x+y)) %>% D("x")## 3 + 1/(x + y)(d)
quote(x*exp(-2*x) + log(1/x)) %>% D("x")## exp(-2 * x) - x * (exp(-2 * x) * 2) - 1/x^2/(1/x)Questão 15
(a)
quote(3*x^3 - cos(x)) %>% D("x") %>% D("x")## 3 * (3 * (2 * x)) + cos(x)(b)
quote(3*x^3 - y*cos(x)) %>% D("x") %>% D("y")## sin(x)Lista 03
library(ggplot2)
library(gridExtra)
library(dplyr)Questão 1
(a)
unif <- sample(0:10, size = 100, replace=T)(b)
gam <- rgamma(100, shape = 2, rate = 3)(c)
norm <- rnorm(100, mean = 3, sd = 9)(d)
expo <- rexp(100, rate = 7)Questão 2
grafb <- ggplot(data = data.frame(x=gam), aes(x=x)) +
geom_histogram(aes(y=..density..), bins = 10) +
stat_function(fun = dgamma, args = list(shape = 2, rate = 3), color = "red") +
theme_minimal() +
labs(title = "Distribuição Gama com shape=2 e rate=3", y="f(x)")
grafc <- ggplot(data = data.frame(x=norm), aes(x=x)) +
geom_histogram(aes(y=..density..), bins = 10) +
stat_function(fun = dnorm, args = list(mean = 3, sd = 9), color = "red") +
theme_minimal() +
labs(title = "Distribuição Normal com media=3 e desv. padrao=9", y="f(x)")
grafd <- ggplot(data = data.frame(x=expo), aes(x=x)) +
geom_histogram(aes(y=..density..), bins = 10) +
stat_function(fun = dexp, args = list(rate = 7), color = "red") +
theme_minimal() +
labs(title = "Distribuição Exponencial com rate=7", y="f(x)")
grid.arrange(grafb, grafc, grafd, ncol=2)
Questão 3
(a)
ppois(0, lambda = 5, lower.tail = F)## [1] 0.9932621(b)
ppois(11, lambda = 5) - ppois(2, lambda = 5)## [1] 0.8698949(c)
qpois(0.76, lambda = 5)## [1] 6(d)
qpois(0.73, lambda = 5, lower.tail = F)## [1] 4## ou
qpois(0.27, lambda = 5)## [1] 4Questão 4
(a)
pbeta(0.7, 3, 10, lower.tail = F)## [1] 0.0002063763(b)
pbeta(0.25, 3, 10) - pbeta(0.2, 3, 10)## [1] 0.1676707(c)
qbeta(0.20, 3, 10)## [1] 0.1307157(d)
qbeta(0.12, 3, 10, lower.tail = F)## [1] 0.3703205## ou
qbeta(0.88, 3, 10)## [1] 0.3703205Questão 5
base <- data.frame(
am.norm = rnorm(200, mean = 10, sd = 5),
am.gama = rgamma(200, rate = 2, shape = 4),
am.pois = rpois(200, lambda = 5),
am.qui = rchisq(200, df = 2)
)(a)
ggplot(base, aes(sample=am.norm)) +
stat_qq(distribution = qchisq, dparams = list(df=2)) +
stat_qq_line(distribution = qchisq, dparams = list(df=2))
qqplot(
x = rchisq(200, 2),
y = rnorm(200, mean = 10, sd = 5)
)
(b)
ggplot(base, aes(sample=am.gama)) +
stat_qq(distribution = qchisq, dparams = list(df=2)) +
stat_qq_line(distribution = qchisq, dparams = list(df=2))
qqplot(
x = rchisq(200, 2),
y = rgamma(200, 3, 4)
)
(c)
ggplot(base, aes(sample=am.pois)) +
stat_qq(distribution = qchisq, dparams = list(df=2)) +
stat_qq_line(distribution = qchisq, dparams = list(df=2))
qqplot(
x = rchisq(200, 2),
y = rpois(200, 2)
)
(d)
ggplot(base, aes(sample=am.qui)) +
stat_qq(distribution = qchisq, dparams = list(df=2)) +
stat_qq_line(distribution = qchisq, dparams = list(df=2))
qqplot(
x = rchisq(200, 2),
y = rchisq(200, 2)
)
Questão 6
(a)
pexp(5, 5, lower.tail = F)## [1] 1.388794e-11(b)
pexp(9, 5) - pexp(2, 5)## [1] 4.539993e-05(c)
dexp(0, 5)## [1] 5(d)
pexp(2, 5)## [1] 0.9999546(e)
dbinom(2, 15, 0.3)## [1] 0.09156011## ou
pbinom(2, 15, 0.3) - pbinom(1, 15, 0.3)## [1] 0.09156011(f)
pbinom(4, 15, 0.3)## [1] 0.5154911(g)
pbinom(5, 15, 0.3) - pbinom(-2, 15, 0.3)## [1] 0.7216214(h)
pbinom(9, 15, 0.3) - pbinom(2, 15, 0.3)## [1] 0.8695198(i)
pbinom(3, 15, 0.3, lower.tail = F)## [1] 0.7031321## ou
1 - pbinom(3, 15, 0.3)## [1] 0.7031321Questão 7
fun7 <- function(y) {
if(y<1) return(0)
if(y>3) return(1)
return(integrate(
(function(x) (x^3)/20), lower = 1, upper = y
)[[1]]
)
}(a)
1 - fun7(2)## [1] 0.8125(b)
fun7(2.6) - fun7(2)## [1] 0.37122(c)
1 - fun7(2.8)## [1] 0.24418(d)
fun7(2) + (1 - fun7(2.5))## [1] 0.7117188Questão 8
(a)
ggplot(data = data.frame(x=c(-1,11)), aes(x=x)) +
stat_function(fun = punif, args = list(min = 0, max = 10)) +
labs(y="F(x)", title="Distribuição acumulada de uma Uniforme(0,10)") +
theme_minimal()
(b)
ggplot(data = data.frame(x=c(0,5)), aes(x=x)) +
stat_function(fun = pgamma, args = list(shape = 2 , rate = 3)) +
labs(y="F(x)", title="Distribuição acumulada de uma Gama(2,3)") +
theme_minimal()
(c)
ggplot(data = data.frame(x=c(-30,40)), aes(x=x)) +
stat_function(fun = pnorm, args = list(mean = 3, sd = 9)) +
labs(y="F(x)", title="Distribuição acumulada de uma N(3,9)") +
theme_minimal()
(d)
ggplot(data = data.frame(x=c(0,2)), aes(x=x)) +
stat_function(fun = pexp, args = list(7)) +
labs(y="F(x)", title="Distribuição acumulada de uma Exponecial(7)") +
theme_minimal()
(e)
ggplot(data = data.frame(x=c(0,4)), aes(x=x)) +
stat_function(fun = Vectorize(fun7)) +
labs(y="F(x)", title="Distribuição acumulada de x^3/20") +
theme_minimal()
(f)
Professor desconsiderou a questão
Questão 9
(a)
ggplot(data = data.frame(x=c(-1:110)), aes(x=x)) +
# geom_point(aes(y=dbinom(x, 100, 0.85)), color="red") +
stat_function(fun = dbinom, args = list(100, 0.85), geom = "point", n = 112)
## ou
plot(dbinom(-1:110,100,0.85), type = 'h')
ggplot(data = data.frame(x=c(-1:110)), aes(x=x)) +
stat_function(fun = pbinom, args = list(100, 0.85), geom = "step") +
labs(y="F(x)")
## ou
plot(pbinom(-1:110,100,0.85), type = 's')
Professor desconsiderou a questão
(b)
Professor desconsiderou a questão
(c)
Professor desconsiderou a questão
Questão 10
Professor desconsiderou a questão
Questão 11
Professor desconsiderou a questão
Questão 12
fun12 <- function(x, k) {
return((gamma(k)*exp(-k*x-3))/k * I(x>0 & k>0))
}(a)
ggplot(data = data.frame(x=c(-2,6)), aes(x=x)) +
stat_function(fun = fun12, args = list(k=4)) +
labs(y="f(x)") +
theme_minimal()
(b)
ggplot(data = data.frame(k=c(-2,6)), aes(x=k)) +
stat_function(fun = fun12, args = list(x=2)) +
labs(y="f(k)") +
theme_minimal()
(c)
maximo <- optimise(fun12, c(-2,6), k=2, maximum = TRUE); maximo## $maximum
## [1] 5.786925e-05
##
## $objective
## [1] 0.02489065(d)
ggplot(data = data.frame(x=c(-2,6)), aes(x=x)) +
stat_function(fun = fun12, args = list(k=2)) +
labs(y="f(x)") +
geom_vline(aes(xintercept = maximo$maximum), linetype = "dashed", color = "red") +
theme_minimal()
Questão 13
(a)
x = c(1, 2, 2.4, 2.8, 5, 1, 3, 4, 6.3, 2.9)
mv_a <- function(beta, x) {
produto = 1
for (i in 1:length(x)) {
produto = produto*(x[i]^2)
}
b <- beta^(3*length(x))
g <- gamma(3)^length(x)
return(((b/g) * exp(-beta*sum(x)) * produto))
}
EMV <- optimise(mv_a, c(0,2), x = c(1, 2, 2.4, 2.8, 5, 1, 3, 4, 6.3, 2.9), maximum = T)
ggplot(data=data.frame(x=c(0,2)), aes(x=x)) +
stat_function(fun = mv_a, args = list(x = c(1, 2, 2.4, 2.8, 5, 1, 3, 4, 6.3, 2.9))) +
geom_vline(xintercept = EMV$maximum, linetype = "dashed", color = "red") +
theme_minimal()
(b)
x = c(10, 12, 11, 12.8, 13, 14.9, 12, 16.2)
mv_b <- function(u,x) {
soma=0
for (i in 1:length(x)) {
soma = soma + ((x[i]-u)/sqrt(10))^2
}
return(((1/sqrt(20*pi))^(length(x))) * exp((-1/2) * soma))
}
EMV <- optimise(mv_b, c(5,20), x = c(10, 12, 11, 12.8, 13, 14.9, 12, 16.2), maximum = T)
ggplot(data = data.frame(x=c(5,20)), aes(x=x)) +
stat_function(fun = mv_b, args = list(x = c(10, 12, 11, 12.8, 13, 14.9, 12, 16.2))) +
geom_vline(xintercept = EMV$maximum, linetype = "dashed", color = "red") +
labs(y="") +
theme_minimal()
(c)
Está indo pro infinito do infinito ..
# x = c(2, 3, 2, 2, 3, 3, 3, 4, 4, 5, 7, 2, 2, 2)
#
# mv_c <- function(p, x) {
# return(((1-p)^sum(x)) * p^length(x))
# }
#
# EMV <- optimise(mv_c, c(300,470), x = c(2, 3, 2, 2, 3, 3, 3, 4, 4, 5, 7, 2, 2, 2), maximum = T)
#
# ggplot(data = data.frame(x=c(300,470)), aes(x=x)) +
# stat_function(fun = mv_c, args = list(x = c(2, 3, 2, 2, 3, 3, 3, 4, 4, 5, 7, 2, 2, 2))) +
# geom_vline(xintercept = EMV$maximum, args = list(x = c(2, 3, 2, 2, 3, 3, 3, 4, 4, 5, 7, 2, 2, 2)), linetype = "dashed", color = "red") +
# theme_minimal()(d)
x = c(25, 23, 22, 21, 27, 39, 35, 33, 32)
mv_d <- function(o2, x) {
soma=0
for (i in 1:length(x)) {
soma = soma + ((x[i]-10)/sqrt(o2))^2
}
return(((1/sqrt(2*pi*o2))^length(x)) * exp((-1/2) * soma))
}
EMV <- optimise(mv_d, c(0,1000), x = c(25, 23, 22, 21, 27, 39, 35, 33, 32), maximum = T)
ggplot(data = data.frame(x=c(0,2000)), aes(x=x)) +
stat_function(fun = mv_d, args = list(x = c(25, 23, 22, 21, 27, 39, 35, 33, 32))) +
geom_vline(xintercept = EMV$maximum, linetype = "dashed", color = "red") +
theme_minimal()
(e)
x <- c(4, 5, 6.2, 4, 3, 5, 6.9, 7, 9.3)
mv_e <- function(lambda, x) {
for (i in 1:length(x)) {
if(x[i]<=0) x[i] <- 0
}
return((lambda^length(x))*exp(-lambda*sum(x)))
}
EMV <- optimise(mv_e, c(0,1), x=c(4, 5, 6.2, 4, 3, 5, 6.9, 7, 9.3), maximum = TRUE)
ggplot(data = data.frame(x=c(0,1)), aes(x=x)) +
stat_function(fun = mv_e, args = list(x = c(4, 5, 6.2, 4, 3, 5, 6.9, 7, 9.3))) +
geom_vline(xintercept = EMV$maximum, linetype = "dashed", color= "red") +
labs(y="") +
theme_minimal() 
Questão 14
P U L A
Questão 15
X <- rnorm(1000,25,sqrt(70))
T <- rchisq(1000, 5)
Y <- rexp(1000,10)
Z <- rgamma(1000,3,5)
W <- rgamma(1000,10,5)(a)
data.frame(x=Z+W) %>%
ggplot(aes(x=x)) +
geom_histogram(aes(y=..density..), bins=30) +
stat_function(fun = dgamma, args = list(13,5), color = "red") +
labs(y="") +
theme_minimal()
(b)
data.frame(x=2*Y) %>%
ggplot(aes(x=x)) +
geom_histogram(aes(y=..density..), bins = 30) +
stat_function(fun = dgamma, args = list(2,10), color = "red") +
labs(y="") +
theme_minimal()
(c)
X_ <- 5*(replicate(1000,rnorm(5,25,sqrt(70)))/5)
mean_X <- apply(X_, 2, mean)
data.frame(x=(mean_X-25)/sqrt(70/5)) %>%
ggplot(aes(x=x)) +
geom_histogram(aes(y=..density..), bins = 30) +
stat_function(fun = dnorm, args = list(0,1), color = "red") +
labs(y="") +
theme_minimal()
(d)
data.frame(x=3*T) %>%
ggplot(aes(x=x)) +
geom_histogram(aes(y=..density..), bins = 30) +
stat_function(fun = dgamma, args = list((15/2), (1/2)), color = "red") +
labs(y="") +
theme_minimal()
(e)
data.frame(x=(X+20) + (X+20) + (X+20)) %>%
ggplot(aes(x=x)) +
geom_histogram(aes(y=..density..), bins = 30) +
stat_function(fun = dnorm, args = list(95, sqrt(210)), color = "red")
Lista 04
library(dplyr)
library(magrittr)
library(readr)
library(ggplot2)Questão 1
(a)
IC_q1a <- function(x, sigma2, confianca) {
amostra = na.omit(x)
n = length(amostra)
x_barra = mean(amostra)
alfa = 1-confianca
z = qnorm(1-alfa)
L = (x_barra - z*sqrt(sigma2/n)) %>% round(1)
cat("IC(m,", confianca*100, "%) = [", L, ",Inf).", sep = "")
}(b)
IC_q1b <- function(x, sigma2, confianca) {
amostra = na.omit(x)
n = length(amostra)
x_barra = mean(amostra)
alfa = 1-confianca
z = qnorm(1-alfa)
L = (x_barra + z*sqrt(sigma2/n)) %>% round(1)
cat("IC(m,", confianca*100, "%) = (-Inf,", L, "].", sep = "")
}(c)
IC_q1c <- function(x, sigma, confianca) {
amostra = na.omit(x)
n = length(amostra)
x_barra = mean(amostra)
alfa = 1-confianca
z = qnorm(1-alfa/2)
LI = (x_barra - z*sigma/sqrt(n)) %>% round(1)
LS = (x_barra + z*sigma/sqrt(n)) %>% round(1)
cat("IC(m,", confianca*100, "%) = [", LI, ",", LS, "].", sep = "")
}(d)
IC_q1d <- function(x, sigma, confianca) {
amostra = na.omit(x)
n = length(amostra)
x_barra = mean(amostra)
alfa = 1-confianca
z = qnorm(1-alfa/2)
LI = (x_barra - z*sigma/sqrt(n)) %>% round(1)
LS = (x_barra + z*sigma/sqrt(n)) %>% round(1)
# cat("IC(m,", confianca*100, "%) = [", LI, ",", LS, "].", sep = "")
cont = ifelse(30>=LI & 30<=LS, 1, 0)
return(cont)
}
amostra = data.frame(
replicate(1000, rnorm(25, mean = 30, sd = sqrt(25)))
)
soma = 0
for (i in 1:length(amostra)) {
ifelse(
(IC_q1d(amostra[[i]], 5, 0.9)==1), (soma = soma +1), 0
)
}
# soma
cat(
(soma/1000)*100, "%", sep = ""
)## 89.9%## Como pode ser observado, a porcentagem é aproximadamente 90%, o que indica que em 90% das vezes o intervalo contém o verdadeiro valor de muQuestão 2
IC_q2 <- function(x, confianca) {
x = na.omit(x)
n = length(x)
p_chapeu = 1/2
alfa = 1 - confianca
z = qnorm(confianca + alfa/2)
erro = sqrt(0.5*0.5/n)*z
LI = p_chapeu - erro
LS = p_chapeu + erro
cat("IC(p,", confianca*100, "%) = [", LI, ",", LS, "]", sep= "")
}Questão 3
base_3 <- read_table2("Lista 04/Base saude.txt", na = "9")
base_3$HIV %<>% factor(labels = c("Não", "Sim"))
base_3$DST %<>% factor(labels = c("Não", "Sim"))
base_3$Sexo %<>% factor(labels = c("Feminino", "Masculino"))Como a função da letra d era uma adaptação da c, usei a função criada na letra c.
(a)
n = length(base_3$Peso); n## [1] 30## Como n é suficientemente grande, pelo TCL podemos estimar a variância pela variância amostral
IC_q1c(base_3$Peso, sigma = sqrt(var(base_3$Peso, na.rm = T)), confianca = 0.97)## IC(m,97%) = [138.5,371.1].## Foi realizada a suposição de que n é suficientemente grande para usar o TCL.(b)
base_q3b <- base_3 %>% filter(HIV=="Sim") %$% Estatura
## Assumindo que a variância da amostra é um bom estimador para a variância da população
IC_q1c(base_q3b, sigma = sqrt(var(base_q3b, na.rm = T)), confianca = 0.95)## IC(m,95%) = [76.7,139.5].## Isso quer dizer que aproximadamente em 95% das vezes o intervalo conterá o verdadeiro valor de mu.
## Foi realizada a suposição de que a variância amostral é um bom estimador para a variância populacional.(c)
Análogo a b
(d)
base_q3d <- base_3 %>% filter(Sexo=="Feminino", HIV=="Não") %$% Idade
## Assumindo que a variância amostral é um bom estimador para a variância populacional
IC_q1c(base_q3d, sigma = sqrt(var(base_q3d, na.rm = T)), confianca = 0.90)## IC(m,90%) = [26.4,42.8].## Isso quer dizer que em 90% das vezes o intervalo conterá o verdadeiro valor de mu.
## Foi realizado a suposição de que a variância amostral é um bom estimador para a variância populacional.Questão 4
IC_q7 <- function(x, y, conf){
grupos = levels(y)
for(i in 1:length(grupos)){
aux = x[y==grupos[i]] %>% .[!is.na(.)]
cat(grupos[i]," - ", "")
IC_q1c(aux, sigma = sqrt(var(aux, na.rm = T)), conf)
cat("\n")
}
}
IC_q7(base_3$Peso, base_3$Sexo, 0.90)## Feminino - IC(m,90%) = [59.1,312.1].
## Masculino - IC(m,90%) = [183.4,423.9].Questão 5
IC_q7(base_3$Estatura, base_3$DST, conf = 0.95)## Não - IC(m,95%) = [56.2,147.2].
## Sim - IC(m,95%) = [97,158.4].## Para os pacientes observados divididos em grupo possui e não possui DST, pode-se observar o intervalo que conterá o verdadeiro valor da média da estatura de cada grupo com um nível de confiança de 95%.
## Com base nos intervalos dos grupos, leva-nos a supor que a média da Estatura de quem possui DST é maior se comparado com a média de quem não possui DST.Questão 6
IC_q6_var <- function(x, conf) {
x = na.omit(x)
n = length(x)
alfa = 1-conf
S2 = var(x, na.rm = T)
x1 = qchisq(alfa/2, df = n-1, lower.tail = F)
x2 = qchisq(1-alfa/2, df = n-1, lower.tail = F)
LI = (n-1)*S2/x1
LS = (n-1)*S2/x2
cat(
"IC(sigma^2,", conf*100, "%) = [", LI %>% round(2), ",", LS %>% round(2), "].", sep = ""
)
}(a)
base_3 %>% filter(Sexo=="Masculino", HIV=="Sim") %>% select(Idade) %>% pull %>% IC_q6_var(0.90)## IC(sigma^2,90%) = [28.5,102.72].## Se realizado o intervalo de confiança para a variância para todas as possíveis amostras da população, teriamos que aproximadamente 90% desses intervalos conteriam o real valor do parâmetro.
## Foram realizadas as suposições de que a amostra veio de uma distribuição normal e que a variância amostral é um bom estimador para a variância populacional.(b)
IC_q6_var(base_3$Estatura, 0.95)## IC(sigma^2,95%) = [3076.14,8934.48].## Se realizado o intervalo de confiança para a variância para todas as possíveis amostras da população, teríamos que aproximadamente 90% desses intervalos conteriam o real valor do parâmetro.
## Foram realizadas as suposições de que a amostra veio de uma distrivui(c)
base_3 %>% filter(Escol==0)## # A tibble: 1 x 10
## Codigo Datacol Sexo Idade Peso Estatura HIV Escol DST Tipo
## <chr> <chr> <fct> <int> <dbl> <dbl> <fct> <int> <fct> <int>
## 1 AB15 01/01/14 Masculino 18 65 166 Sim 0 Não NA## Não é possível encontrar o intervalo de confiança pois o tamanho da amostra é 1.Questão 7
base_7 <- read_rds("Lista 04/exames medicos.rds")(a)
## HDL
## Empiricamente, não é razoável supor normalidade.
ggplot(base_7, aes(sample = base_7$HDL)) + geom_qq() + geom_qq_line()
## LDL
## Empiricamente, é razoável supor normalidade.
ggplot(base_7, aes(sample = base_7$LDL)) + geom_qq() + geom_qq_line()
## Glicose
## Empiricamente, é razoável supor normalidade
ggplot(base_7, aes(sample = base_7$glicose)) + geom_qq() + geom_qq_line()
## Linfocitos
## Empiricamente, não é razoável supor normalidade.
ggplot(base_7, aes(sample = base_7$linfocitos)) + geom_qq() + geom_qq_line()
(b)
IC_media_q7 <- function(x, conf){
x_barra = mean(x, na.rm = T)
s2 = var(x, na.rm = T)
n = length(x[!is.na(x)])
alfa = 1-conf
erro = qt(1-alfa/2, df = n-1) * sqrt(s2/n)
LI = x_barra - erro
LS = x_barra + erro
cat(
"IC (m,", conf*100, "%) = [", LI, ";", LS, "]", sep = ""
)
}IC_media_q7(base_7$LDL, 0.98)## IC (m,98%) = [29.06364;31.40836]IC_media_q7(base_7$glicose, 0.98)## IC (m,98%) = [46.36625;51.08975](c)
## Empiricamente. Igual feito na letra a.(d)
IC_var_q7 <- function(x, conf){
n = length(x[!is.na(x)])
s2 <- var(x, na.rm = T)
alfa = 1-conf
x1 = qchisq(alfa/2, df = n-1, lower.tail = F)
x2 = qchisq(1-alfa/2, df = n-1, lower.tail = F)
LI = (n-1)*s2 / x1
LS = (n-1)*s2 / x2
cat(
"IC (sigma,", conf*100, "% = [", LI, ";", LS, "]", sep = ""
)
}IC_var_q7(base_7$LDL, 0.95)## IC (sigma,95% = [3.373065;10.70688]IC_var_q7(base_7$glicose, 0.95)## IC (sigma,95% = [13.68889;43.45165](e)
IC_prop_q7 <- function(num_s, n, conf, conservador=FALSE){
if_else(
conservador==FALSE, p_chapeu = num_s/n, p_chapeu = 1/2
)
alfa = 1-conf
erro = qnorm(1-alfa/2) * sqrt( (pchapeu*(1-p_chapeu))/n )
LI = p_chapeu - erro
LS = p_chapeu + erro
cat(
"IC (p,", conf*100, "%) = [", LI, ";", LS, "]", sep = ""
)
}Questão 8
base_8 <- read_rds("Lista 04/colesterol.rds")(a)
## Empiricamente, não é razoável supor normalidade.
ggplot(base_8, aes(sample = base_8$HDL)) + geom_qq() + geom_qq_line()
## Empiricamente, é razoável supor normalidade.
ggplot(base_8, aes(sample = base_8$LDL)) + geom_qq() + geom_qq_line()
(b)
## Suponde normalidade independente do resultado na questão acima.
## Supondo normalidade, sim é possível obter o Intervalo de Confiança para a média.
IC_media_q8 <- function(x, conf){
n = length(x[!is.na(x)])
x_barra = mean(x, na.rm = T)
s2 = var(x, na.rm = T)
alfa = 1-conf
erro = qt(1-alfa, df = n-1) * sqrt(s2/n)
LI = x_barra - erro
LS = x_barra + erro
cat(
"IC (m,", conf*100, ") = [", LI, ";", LS, "]", sep = ""
)
}## Isso quer dizer que aproximadamente 95% das vezes o intervalo conterá o verdadeiro valor de mu.
base_8$HDL %>% IC_media_q8(0.95)## IC (m,95) = [0.2874198;0.3948024]## Isso quer dizer que aproximadamente 95% das vezes o intervalo conterá o verdadeiro valor de mu.
base_8$LDL %>% IC_media_q8(0.95)## IC (m,95) = [49.69716;50.44951](c)
## Empiricamente, igual feito na letra a.(d)
## Suponde normalidade independente do resultado na questão acima.
## Supondo normalidade, sim é possível obter o Intervalo de Confiança para a variância.
IC_var_q8 <- function(x, conf){
n = length(x[!is.na(x)])
s2 = var(x, na.rm = T)
alfa = 1-conf
x1 = qchisq(alfa/2, df = n-1)
x2 = qchisq(1-alfa/2, df = n-1)
LI = (n-1)*s2/x2
LS = (n-1)*s2/x1
cat(
"IC (sigma,", conf*100, ") = [", LI, ";", LS, "]", sep = ""
)
}## Isso quer dizer que aproximadamente 95% das vezes o intervalo conterá o verdadeiro valor de sigma.
base_8$HDL %>% IC_var_q8(0.95)## IC (sigma,95) = [0.07144162;0.128993]## Isso quer dizer que aproximadamente 95% das vezes o intervalo conterá o verdadeiro valor de sigma.
base_8$LDL %>% IC_var_q8(0.95)## IC (sigma,95) = [3.50696;6.332068]Lista 05
library(dplyr)
library(ggplot2)Questão 1
(a)
## Amostras a. s. obtidas de uma população normal.
## Variância conhecida.
## H0: mu = mu0
## H1: mu < mu0
TH_media_q1a <- function(x, mu0, sigma2, alfa){
amostra = na.omit(x)
n = length(amostra)
z = qnorm(1-alfa)
x_barra = mean(amostra)
x_critico = mu0 - z*sqrt(sigma2/n)
decisao = ifelse((x_barra < x_critico), "Rejeita H0.", "Não rejeita H0.")
cat(
"Estimativa pontual: ", x_barra, "\n",
"Região Crítica: [", -Inf, ";", x_critico, "] \n",
"Decisão: ", decisao,
sep = ""
)
}
## mu0 = 500
rnorm(1000, 500, sqrt(25)) %>% TH_media_q1a(500, 25, 0.05)## Estimativa pontual: 500.0001
## Região Crítica: [-Inf;499.7399]
## Decisão: Não rejeita H0.(b)
## Amostras a. s. obtidas de uma população normal.
## Variância desconhecida.
## H0: sigma2 = sigma0
## H1: sigma2 > sigma0
TH_sigma2_q1b <- function(x, sigma0, alfa){
amostra = na.omit(x)
n = length(amostra)
s2 = var(amostra)
x = (n-1)*s2 / sigma0
x_critico = qchisq(alfa, n-1, lower.tail = F)
decisao = ifelse(sigma0 > x_critico, "Rejeita H0.", "Não rejeita H0.")
cat(
"Estimativa pontual: ", x, "\n",
"Região Crítica: [", x_critico, ";", +Inf, "] \n",
"Decisão: ", decisao,
sep = ""
)
}
## sigma0 = 25
rnorm(1000, 500, sqrt(25)) %>% TH_sigma2_q1b(25, 0.05)## Estimativa pontual: 929.1633
## Região Crítica: [1073.643;Inf]
## Decisão: Não rejeita H0.(C)
## Amostras a. s. obtidas de uma população normal.
## H0:
## H1:
TH_prop_q1c <- function(x, p0, alfa){
amostra = na.omit(x)
n = length(amostra)
q = qnorm(1-alfa/2)
k = sqrt( (p0*(1-p0))/n ) * q
LI = p0 - k
LS = p0 + k
decisao = ifelse(( p0<LI | p0>LS ), "Rejeita H0.", "Não rejeita H0.")
cat(
"Estimativa pontual: ", mean(amostra), "\n",
"Região Crítica: [", LI, ";", LS, "] \n",
"Decisão: ", decisao,
sep = ""
)
}
c(1,1,1,1,1,1,0,0,0,0,1,0,0,1) %>% TH_prop_q1c(., sum(./length(.)), 0.05)## Estimativa pontual: 0.5714286
## Região Crítica: [0.3122037;0.8306534]
## Decisão: Não rejeita H0.(d)
## H0: mu = mu0
## H1: mu < mu0
TH_media_q1d <- function(x, mu0, sigma2, alfa){
amostra = na.omit(x)
n = length(amostra)
x_barra = mean(amostra)
z = (x_barra-mu0) / sqrt(sigma2/n)
q = qnorm(alfa)
decisao = ifelse((z < q), "Rejeita H0.", "Não rejeita H0.")
cat(
"Estimativa pontual: ", x_barra, "\n",
"Estimativa de teste: ", z, "\n",
"Região Crítica: [", -Inf, ";", q, "] \n",
"Decisão: ", decisao,
sep = ""
)
}
rnorm(1000, 500, sqrt(25)) %>% TH_media_q1d(500, 25, 0.05)## Estimativa pontual: 499.9604
## Estimativa de teste: -0.2503271
## Região Crítica: [-Inf;-1.644854]
## Decisão: Não rejeita H0.(e)
## Erro tipo II: Probabilidade de não rej. H0 dado H0 verdadeiro.
erroII_q1e <- function(n, alfa, sigma2, mu0, mu){
# Região Critica
x_critico = mu0 - sqrt(sigma2/n)*qnorm(1-alfa)
# Padronizando
z = (x_critico - mu) / sqrt(sigma2/n)
return( pnorm(z, lower.tail = F) )
}
erroII_q1e(1000, 0.05, 25, 500, 500)## [1] 0.95(f)
## O gráfico mostra as probabilidades de se não rejeitar H0 para cada mu diferente.
ggplot(data = data.frame(x = c(0, 60)), aes(x = x)) +
stat_function(fun = erroII_q1e, args = list(n = 20, alfa = 0.05, sigma2 = 25, mu0 = 30)) +
geom_vline(xintercept = 30, linetype = "dashed", col = "red") +
ggtitle("Erro tipo II") + xlab(expression(mu)) + ylab(expression(beta(mu))) 
(g)
## A medida que o tamanho da amostra aumenta os valores a esquerda de mu0 ficam menos prováveis, fazendo a curva crescer quando está mais próximo de mu0
ggplot(data = data.frame(x = c(0, 60)), aes(x = x)) +
stat_function(fun = erroII_q1e, args = list(n = 20 , alfa = 0.05, sigma2 = 25, mu0 = 30), color = "red") +
stat_function(fun = erroII_q1e, args = list(n = 30 , alfa = 0.05, sigma2 = 25, mu0 = 30), color = "black") +
stat_function(fun = erroII_q1e, args = list(n = 40 , alfa = 0.05, sigma2 = 25, mu0 = 30), color = "purple") +
stat_function(fun = erroII_q1e, args = list(n = 50 , alfa = 0.05, sigma2 = 25, mu0 = 30), color = "green") +
stat_function(fun = erroII_q1e, args = list(n = 100, alfa = 0.05, sigma2 = 25, mu0 = 30), color = "blue") +
geom_vline(xintercept = 30, linetype = "dashed", col ="red") +
ggtitle("Erro tipo II") + xlab(expression(mu)) + ylab(expression(beta(mu))) 
(h)
poder_q1h <- function(n, alfa, sigma2, mu0, mu){
# Região Crítica
q = qnorm(1-alfa/2)
LI = mu0 - sqrt(sigma2/n)*q
LS = mu0 + sqrt(sigma2/n)*q
# Padronizando
z1 = (LI - mu) / sqrt(sigma2/n)
z2 = (LS - mu) / sqrt(sigma2/n)
return(pnorm(z1) + pnorm(z2, lower.tail = F))
}
rnorm(1000, 500, sqrt(25)) %>% length() %>% poder_q1h(0.05, 25, 500, 500)## [1] 0.05(i)
## O gráfico mostra as probabilidades de se rejeitar H0 para cada mu diferente.
ggplot(data = data.frame(x = c(0, 60)), aes(x = x)) +
stat_function(fun = poder_q1h, args = list(n = 20, alfa = 0.05, sigma2 = 25, mu0 = 30)) +
geom_vline(xintercept = 30, linetype = "dashed", col ="red") +
ggtitle("Função Poder") + xlab(expression(mu)) + ylab(expression(pi(mu))) 
(j)
## A medida que o tamanho da amostra aumenta maior é o poder do teste, ou seja, a probabilidade de se rejeitar valores diferentes de mu0 aumenta.
ggplot(data = data.frame(x = c(0, 60)), aes(x = x)) +
stat_function(fun = poder_q1h, args = list(n = 20 , alfa = 0.05, sigma2 = 25, mu0 = 30), color = "red") +
stat_function(fun = poder_q1h, args = list(n = 30 , alfa = 0.05, sigma2 = 25, mu0 = 30), color = "black") +
stat_function(fun = poder_q1h, args = list(n = 40 , alfa = 0.05, sigma2 = 25, mu0 = 30), color = "purple") +
stat_function(fun = poder_q1h, args = list(n = 50 , alfa = 0.05, sigma2 = 25, mu0 = 30), color = "green") +
stat_function(fun = poder_q1h, args = list(n = 100, alfa = 0.05, sigma2 = 25, mu0 = 30), color = "blue") +
geom_vline(xintercept = 30, linetype = "dashed", col ="red") +
ggtitle("Função Poder") + xlab(expression(mu)) + ylab(expression(pi(mu))) 
Questão 2
(a)
## Amostras a. s. obtidas de uma população normal com média mu e variância sigma2 desconhecida.
## Suponha que n é grande.
## H0: mu = mu0
## H1: mu > mu0
TH_media_q2a <- function(x, alfa, mu0){
amostra = na.omit(x)
n = length(amostra)
## x_barra segue distribuição normal pois é combinação linear de distribuições normais.
x_barra = mean(amostra)
s2 = var(amostra)
q = qnorm(alfa, lower.tail = F)
k = sqrt(s2/n) * q
x_critico = mu0 + k
decisao = if_else((x_barra > x_critico), "Rejeita H0", "Não rejeita H0")
cat(
"Estimativa pontual: ", x_barra, "\n",
"Região Crítica: [", x_critico, " ; ", +Inf, "] \n",
"Decisão: ", decisao, " a um nível de significância ", alfa,
sep = ""
)
}
rnorm(1000, 500, sqrt(25)) %>% TH_media_q2a(0.05, 500)## Estimativa pontual: 499.8381
## Região Crítica: [500.2596 ; Inf]
## Decisão: Não rejeita H0 a um nível de significância 0.05(b)
## Amostras a. s. obtidas de uma população com distribuição normal com média mu e variância desconhecida.
## H0: mu = mu0
## H1: mu > mu0
TH_media_q2b <- function(x, alfa, mu0){
amostra = na.omit(x)
n = length(amostra)
## x_barra segue distribuição normal pois é combinação linear de distribuições normais.
x_barra = mean(amostra)
s2 = var(amostra)
t = (x_barra - mu0) / sqrt(s2/n)
x_critico = qnorm(alfa, lower.tail = F)
decisao = if_else((t > x_critico), "Rejeita H0", "Não rejeita H0")
cat(
"Estatística de teste: ", t, "\n",
"Região Crítica: [", x_critico, " ; ", +Inf, "] \n",
"Decisão: ", decisao, " a um nível de significância ", alfa,
sep = ""
)
}
rnorm(1000, 500, sqrt(25)) %>% TH_media_q2a(0.05, 500)## Estimativa pontual: 500.2115
## Região Crítica: [500.2586 ; Inf]
## Decisão: Não rejeita H0 a um nível de significância 0.05(c)
erro_tipoII_q2c <- function(n, alfa, s2, mu0, mu){
## Região critica
k = sqrt(s2/n)*qnorm(1-alfa)
x_critico = mu0 + k
## Padronizando e calculando probabilidade
z = (x_critico - mu) / sqrt(s2/n)
return( pnorm(z) )
}
erro_tipoII_q2c(1000, 0.05, 25, 500, 500)## [1] 0.95(d)
## O gráfico mostra as probabilidades de se não rejeitar H0 para cada mu diferente.
ggplot(data = data.frame(x = c(0, 20)), aes(x = x)) +
stat_function(fun = erro_tipoII_q2c, args = list(n = 20, alfa = 0.05, s2 = 25, mu0 = 10)) +
geom_vline(xintercept = 10, linetype = "dashed", col = "red") +
ggtitle("Erro tipo II") + xlab(expression(mu)) + ylab(expression(beta(mu)))
(e)
## Pode-se afirmar que a medida que o tamanho da amostra aumenta os valores a direita de mu0 ficam menos prováveis.
ggplot(data = data.frame(x = c(0, 60)), aes(x = x)) +
stat_function(fun = erro_tipoII_q2c, args = list(n = 20 , alfa = 0.05, s2 = 25, mu0 = 30), color = "red") +
stat_function(fun = erro_tipoII_q2c, args = list(n = 30 , alfa = 0.05, s2 = 25, mu0 = 30), color = "black") +
stat_function(fun = erro_tipoII_q2c, args = list(n = 50 , alfa = 0.05, s2 = 25, mu0 = 30), color = "green") +
stat_function(fun = erro_tipoII_q2c, args = list(n = 100, alfa = 0.05, s2 = 25, mu0 = 30), color = "blue") +
geom_vline(xintercept = 30, linetype = "dashed", col = "red") +
ggtitle("Erro tipo II") + xlab(expression(mu)) + ylab(expression(beta(mu)))
(f)
ggplot(data = data.frame(x = c(0, 60)), aes(x = x)) +
stat_function(fun = erro_tipoII_q2c, args = list(n = 20, alfa = 0.01, s2 = 25, mu0 = 30), color = "red") +
stat_function(fun = erro_tipoII_q2c, args = list(n = 20, alfa = 0.05, s2 = 25, mu0 = 30), color = "black") +
stat_function(fun = erro_tipoII_q2c, args = list(n = 20, alfa = 0.07, s2 = 25, mu0 = 30), color = "green") +
stat_function(fun = erro_tipoII_q2c, args = list(n = 20, alfa = 0.10, s2 = 25, mu0 = 30), color = "blue") +
geom_vline(xintercept = 30, linetype = "dashed", col = "red") +
ggtitle("Erro tipo II") + xlab(expression(mu)) + ylab(expression(beta(mu)))
(g)
poder_q2 <- function(n, alfa, s2, mu0, mu){
## Região Crítica
k = sqrt(s2/n)*qnorm(1-alfa)
x_critico = mu0 + k
## Padronizando e calculando probabilidade
z = (x_critico - mu) / sqrt(s2/n)
return(pnorm(z, lower.tail = F))
}
poder_q2(1000, 0.05, 25, 500, 500)## [1] 0.05## Plotando a função
ggplot(data = data.frame(x = c(0,20)), aes(x = x)) +
stat_function(fun = poder_q2, args = list(n = 20, alfa = 0.05, s2 = 25, mu0 = 10)) +
geom_vline(xintercept = 10, linetype = "dashed", col = "red") +
ggtitle("Erro tipo II") + xlab(expression(mu)) + ylab(expression(pi(mu)))
Questão 3
base_3 <- readRDS("Lista 05/BaseGenero.rds")(a)
## H0: mu = 70
## H1: mu > 70
## Sim, com base no nível de significância de 5%, há evidências de que o pesquisador está certo.
base_3 %>% filter(Sexo=="Homem") %>% select(Peso) %>%
t.test(alternativa = "greater", mu = 70) ##
## One Sample t-test
##
## data: .
## t = 13.065, df = 24, p-value = 2.107e-12
## alternative hypothesis: true mean is not equal to 70
## 95 percent confidence interval:
## 79.31101 82.80456
## sample estimates:
## mean of x
## 81.05778(b)
## H0: mu = 30
## H1: mu < 30
## Não, com base no nível de significância de 1% há evidências de que o pesquisador está errado.
TH <- function(x, mu0, alfa){
amostra = na.omit(x)
n = length(amostra)
## Supondo que n é grande o suficiente
s2 = var(amostra)
k = sqrt(s2/n) * qt(alfa, n-1, lower.tail = F)
x_critico = mu0 - k
cat(
"RC = [", -Inf, " ; ", x_critico, "]", sep = ""
)
}
base_3 %>% filter(Sexo=="Homem") %>% select(Idade) %>% pull %>% TH(30, 0.01)## RC = [-Inf ; 27.73446](c)
## H0: mu = 20
## H1: mu =/= 20
## Rejeita h0; Sim, com base no nível de significância de 3% há evidências de que o pesquisador está correto.
base_3 %>% filter(Sexo=="Mulher") %>% select(Idade) %>% t.test(conf.level = 0.99, alternativa = "two.sided", mu = 20)##
## One Sample t-test
##
## data: .
## t = 6.5822, df = 19, p-value = 2.667e-06
## alternative hypothesis: true mean is not equal to 20
## 99 percent confidence interval:
## 22.96809 27.53191
## sample estimates:
## mean of x
## 25.25(d)
## H0: p = 65
## H1: p > 65
## Com base nos dados e adotando um nível de significância de 5%, não rejeitamos H0, ou seja, há evidência de que a proporção de mulheres com peso superior a 70kg não é superior a 65%.
amos <- base_3 %>% filter(Sexo=="Mulher") %>% select(Peso) %>% pull
proporcao = ifelse(amos>70, 1, 0)
prop.test(sum(proporcao), length(proporcao), alternative = "greater", p = 0.65, correct = F)##
## 1-sample proportions test without continuity correction
##
## data: sum(proporcao) out of length(proporcao), null probability 0.65
## X-squared = 14.066, df = 1, p-value = 0.9999
## alternative hypothesis: true p is greater than 0.65
## 95 percent confidence interval:
## 0.1273772 1.0000000
## sample estimates:
## p
## 0.25Questão 4
(a)
t.test()
(b)
t.test()
(c)
t.test()
(d)
prop.test()
Questão 5
erro_tipoII <- function(n, p0, p, alfa){
## Região Crítica
pzin = sqrt(p0*(1-p0)/n) * qnorm(alfa, lower.tail = F)
p_critico = p0 + pzin
## Probabilidade
return(
1-pnorm(p_critico, p, sqrt(p*(1-p)/n), lower.tail = F )
)
}## Gráfico da função
ggplot(data.frame(x = c(0,1)), aes(x = x)) +
stat_function(fun = erro_tipoII, args = list(n = 20, p0 = 0.65, alfa = 0.05))
Questão 6
poder <- function(n, p0, p, alfa){
## Região Crítica
pzin = sqrt(p0*(1-p0)/n) * qnorm(alfa, lower.tail = F)
LI = p0 - pzin
LS = p0 + pzin
## Probabilidade
z1 = (LI - p)/sqrt(p*(1-p)/n)
z2 = (LS - p)/sqrt(p*(1-p)/n)
return(
pnorm(z1) + pnorm(z2, lower.tail = F)
)
}## Gráfico da função
ggplot(data.frame(x = c(0,1)), aes(x = x)) +
stat_function(fun = poder, args = list(n = 20, p0 = 0.65, alfa = 0.05))