Prova 2 de EC
Edimar Fernandes
2024-04-10
library(DT)
library(dplyr)
library(ggplot2)
library(DAAG)1.Considerando o cálculo da integral pelo método de Monte Carlo:
\(\int_{0}^{\infty} \log(1 + x)e^{−x}\,dx\)
***a.Simule 10.000 realizações da variável aleatória exponencial de parâmetro 1 e utilize essas realizações para:
Estimar o valor da integral Obtenha também a estimativa do erro padrão do estimador da integral.
\[ \Large\int_{0}^{\infty} \log(1 + x)e^{-x}\, dx \]
n <-10000
Int.exp <- function(n){
x = rexp(n) # usamos rexp para gerar números aleatórios seguindo uma distribuição exponencial
y <- log(1 + x)*exp(-x) # a formula da função de distruibuição Exponencial
I <- mean(y) # calcula a estimativa do valor da integral usando o método de Monte Carlo.
erro_pad <- sd(y)/ sqrt(n) # Erro padrão
return(list(integral = I, erro = erro_pad)) # criando um objeto para estimativa do valor da integral e o erro padrão do estimador da integral
}
resultado <- Int.exp(n)
print(resultado$integral)## [1] 0.1805614print(resultado$erro)## [1] 0.00076935482.Um dado experimento aleatório consiste em lançar dois dados não viesados (6 faces em cada dado) e observar a soma obtida. Por meio de um procedimento de MC, obtenha a probabilidade aproximada da soma ser par. Simule para N = 10, 100, 1000, 10000 e 100000, em que N é o número de réplicas de MC. O que você observa? Explique.
Criando o procedimento de Monte Carlo
Fazendo simulações com N=10
N = 10
contar_par10 <- 0
dado1 <- 1:6
dado2 <- 1:6
lance1 <- numeric(N)
lance2 <- numeric(N)
for(i in 1:N){
lance1 [i] <- sample(dado1, 1)
lance2 [i] <- sample(dado2, 1)
soma <- lance1[i] + lance2[i]
if (soma %% 2 == 0){
contar_par10 <- contar_par10 + 1
}
}Vamos calcular a probabilidade aproximada da soma ser par.
prob_contar_par10 <- contar_par10/N
cat(" A probabilidade de obter numeros pares nesse experimento é :", prob_contar_par10, "\n")## A probabilidade de obter numeros pares nesse experimento é : 0.6Fazendo simulações com N=100
N = 100
contar_par100<- 0
dado1 <- 1:6
dado2 <- 1:6
lance1 <- numeric(N)
lance2 <- numeric(N)
for(i in 1:N){
lance1 [i] <- sample(dado1, 1)
lance2 [i] <- sample(dado2, 1)
soma <- lance1[i] + lance2[i]
if (soma %% 2 == 0){
contar_par100 <- contar_par100 + 1
}
}
prob_contar_par100 <- contar_par100/N
cat(" A probabilidade de obter numeros pares nesse experimento é :", prob_contar_par100, "\n")## A probabilidade de obter numeros pares nesse experimento é : 0.54Fazendo simulações com N=1000
N = 1000
contar_par1000 <- 0
dado1 <- 1:6
dado2 <- 1:6
lance1 <- numeric(N)
lance2 <- numeric(N)
for(i in 1:N){
lance1 [i] <- sample(dado1, 1)
lance2 [i] <- sample(dado2, 1)
soma <- lance1[i] + lance2[i]
if (soma %% 2 == 0){
contar_par1000 <- contar_par1000 + 1
}
}
prob_contar_par1000 <- contar_par1000/N
cat(" A probabilidade de obter numeros pares nesse experimento é :", prob_contar_par1000, "\n")## A probabilidade de obter numeros pares nesse experimento é : 0.503Fazendo simulações com N=10000
N = 10000
contar_par10000 <- 0
dado1 <- 1:6
dado2 <- 1:6
lance1 <- numeric(N)
lance2 <- numeric(N)
for(i in 1:N){
lance1 [i] <- sample(dado1, 1)
lance2 [i] <- sample(dado2, 1)
soma <- lance1[i] + lance2[i]
if (soma %% 2 == 0){
contar_par10000 <- contar_par10000
}
}
prob_contar_par10000 <- contar_par10000/N
cat(" A probabilidade de obter numeros pares nesse experimento é :", prob_contar_par10000, "\n")## A probabilidade de obter numeros pares nesse experimento é : 0Fazendo simulações com N=100000
N = 100000
contar_par100000 <- 0
dado1 <- 1:6
dado2 <- 1:6
lance1 <- numeric(N)
lance2 <- numeric(N)
for(i in 1:N){
lance1 [i] <- sample(dado1, 1)
lance2 [i] <- sample(dado2, 1)
soma <- lance1[i] + lance2[i]
if (soma %% 2 == 0){
contar_par100000 <- contar_par100000 + 1
}
}
prob_contar_par100000 <- contar_par100000/N
cat(" A probabilidade de obter numeros pares nesse experimento é :", prob_contar_par100000, "\n")## A probabilidade de obter numeros pares nesse experimento é : 0.50201——————————————–Resumo das amostras N—————————————————–
# Criar um data frame com os resultados
amostra_N <- data.frame(
N = c(10, 100, 1000, 10000, 100000),
Probabilidade = c(prob_contar_par10, prob_contar_par100, prob_contar_par1000, prob_contar_par10000, prob_contar_par100000)
)
# Criando uma tabela para visualização do dataframe
datatable(amostra_N,extensions = "Buttons",
options = list(dom = "Bfrtip",
buttons = c("pdf"), pageLength = 5))# Plotando um grafico para ver a evolução:
(graf_Npar <- ggplot(amostra_N, aes(x = N, y = Probabilidade))+
geom_point()+ # adiciona pontos ao gráfico
geom_line()+
labs(x = "Tamanho das amostras calculadas (N)", y = "Probabilidades")+
ggtitle("Simulação para N = 10, 100, 1000, 10000 e 100000, em que N é o número de réplicas de MC")+
xlim(0, max(amostra_N$N)) + # Define os limites do eixo x
ylim(0,1)) # Define os limites do eixo y
# PODEMOS PERCEBER QUE A MEDIDA QUE NOSSA AMOSTRA N AUMENTA TEREMOS UMA ESTIMATIVA MAIS PRECISA DA PROBABILIDADE, OU SEJA , A VARIABILIDADE DE ESTIMATIVA DIMINUI A MEDIDA QUE AUMENTA ESSSA AMOSTRA N!!!3. Considere um dado VICIADO de 6 faces, com faces iguais a: 1, 1, 1, 2, 2 e 3. Ou seja, três faces iguais a 1, duas faces iguais a 2 e uma face igual a 3. Em dois lançamentos independentes desse dado, obtenha a probabilidade de:
i. obter a soma das faces igual a um número par.
ii. a soma ser menor ou igual a 3.
iii.a soma não ser divisível por 3.
iv. Qual é o valor esperado da soma das faces dos dois dados? Apresente a rotina computacional de todos os itens.
# PODEMOS PENSAR EM SEPARAR NOSSOS DADOS EM dois OBJETOS pois são dois lançamentos
n <- 1000 # Repetições de monte Carlo cada vez que aumentamos essa amostra temos a probabilidade mais real proxima
dado <- c(1,1,1,2,2,3) # faces iguais a: 1, 1, 1, 2, 2 e 3
sorteio1 <- numeric(n)
sorteio2 <- numeric(n)
for(i in 1:n){
sorteio1[i] <- sample(dado,1, replace = T)
sorteio2[i] <- sample(dado,1, replace = T)
}
# Mostrando os 10 primeiros sorteios de cada objeto
head(sorteio1)## [1] 2 1 2 1 1 1head(sorteio2)## [1] 2 1 1 3 2 1Resolvendo as probabilidades
# i. Probabilidade de obter a soma das faces igual a um número par
prob_par <- sum((sorteio1 + sorteio2) %% 2 == 0) / n # somando as faces dos dois dados, %% resto da divisão por 2, == 0 verifica se o resto da divisão é igual a zero, se a soma é divisível por 2 significando que é um numero par, sum((sorteio1 + sorteio2) %% 2 == 0) quantas vezes a soma das faces dos dois dados é um número par
cat("Logo temos que a probabilidade de obter a soma das faces igual a um número par:", prob_par, "\n")## Logo temos que a probabilidade de obter a soma das faces igual a um número par: 0.552# ichi. Probabilidade da soma ser menor ou igual a 3
prob_menor_igual_3 <- sum(sorteio1 + sorteio2 <= 3) / n
cat("Logo temos que a probabilidade da soma ser menor ou igual a 3:", prob_menor_igual_3, "\n")## Logo temos que a probabilidade da soma ser menor ou igual a 3: 0.583# iii. Probabilidade da soma não ser divisível por 3
prob_nao_divisivel_3 <- sum((sorteio1 + sorteio2) %% 3 != 0) / n
cat("Logo temos que a probabilidade da soma não ser divisível por 3:", prob_nao_divisivel_3, "\n")## Logo temos que a probabilidade da soma não ser divisível por 3: 0.641# iv. Valor esperado da soma das faces dos dois dados
esperanca <- mean(sorteio1 + sorteio2)
cat("Valor esperado da soma das faces dos dois dados:", esperanca, "\n")## Valor esperado da soma das faces dos dois dados: 3.3344. Exercise 6.1.(p. 178) - Compute a Monte Carlo estimate of
\[ \Large\int_{0}^{\frac{\pi}{3}}\sin(t) \, dt \]
#set.seed(123)
n <- 10000
min <- 0
max <- pi/3
integral <- function(n, min, max){
x <- runif(n, min, max)
prob <- sin(x)
integral_interv <- (max - min)*mean(prob)
return(integral_interv)
}
integral <- integral(n, min, max)
print(integral)## [1] 0.497818Podemos fazer também uma maneira mais direta usando a funçao integrate
set.seed(123)
n = 10000
# Defina a função a ser integrada
I <- function(I, n) {
sin(I)
}
# Intervalo de integração limites inferior e superior de integração.
min <- 0
max <- pi/3
# Calcule a integral
result_integral <- integrate(I, min, max) # Usando a função integrate para calcular a integral da função I
# Exiba o resultado
print(result_integral)## 0.5 with absolute error < 5.6e-155. Explique o Example 7.7. (p. 193).
Example 7.7 (taxa de erro empírico tipo I).
Suponha que \(X_1, . . . ,X_{20}\) é uma amostra aleatória de uma distribuição \(N(\mu, \sigma^2)\).
Teste \(H_0 : \mu = 500_{H_1} : \mu > 500\) em \(\alpha=0,05\) Sob a hipótese nula, \(\large\ T^* =\frac{X- 500}{S\sqrt{20}} \sin{}t(19)\)
onde t(19) denota a distribuição t de Student com 19 graus de liberdade. Grandes valores de \(T^*\) apoiam a hipótese alternativa.
Use um Monte Carlo método para calcular uma probabilidade empírica de erro Tipo I quando $= 100 $, e verificar se é aproximadamente igual a \(\alpha= 0,05.\).
A simulação abaixo ilustra o procedimento para o caso \(\sigma= 100.\). O teste t é implementado por t.test em R, e estamos baseando as decisões de teste nos valores p relatados retornados por t.test.
n <- 20 # amostra aleatória de valores da distribuição NORMAL
alpha <- .05 # nível de significancia
mu0 <- 500 # valor da media da hipótese nula
sigma <- 100 # valor do desvio padrão como 100
m <- 10000 # gerando m realizações de de Conte Carlo
p <- numeric(m) # Cria um vetor p VAZIO e armazenando p valores para cada simulação (1000 vezes)
for (j in 1:m) {# varrendo j de 1 até m, ou seja um loop
x <- rnorm(n, mu0, sigma) # gera n amostras aleatórias (n=20 amostras), a partir da distribuição NORMAL, com média mi e desvio padrão sigma
ttest <- t.test(x, alternative = "greater", mu = mu0) # Realiza um teste unilateral a direita , utiliza a a mostra aletória x e a média mu0 = 500
p[j] <- ttest$p.value # Armazena o valor do teste t
}
p.hat <- mean(p < alpha) # Calcula a probabilidade empírica de erro tipo I. Numero vezes em que o valor p observado é menor que o nível de significância alpha = 5%
se.hat <- sqrt(p.hat * (1 - p.hat) / m) # Calcula o erro padrão da prob Impítica de erro tipo I
print(c(p.hat, se.hat)) # Mostrando os valores de probabilidade empírica de erro tipo e o erro padrão## [1] 0.053800000 0.002256226Vamos criar um cenário das 10000 amostras
# Geração de 10.000 amostras sob a hipótese nula
set.seed(123) # Define uma semente para reprodução dos resultados
amostras <- replicate(10000, { #
x <- rnorm(n, mu0, sigma)
mean(x)
})
# Média da amostra observada
media_amostra <- mean(amostras)
# Limite crítico do teste unilateral à direita
limite_critico <- qnorm(1 - alpha)Vamos fazer um grafico para melhor entendimento visual
# Plotagem do histograma das médias amostrais
hist(amostras, breaks = 30, freq = FALSE, xlab = "Médias Amostrais", ylab = "Densidade",
main = "Distribuição das Médias Amostrais")
abline(v = media_amostra, col = "red", lwd = 2, lty = 2) # Adiciona uma linha vertical para a média da amostra
abline(v = mu0 + limite_critico * sigma / sqrt(n), col = "blue", lwd = 2) # Adiciona uma linha vertical para o limite crítico
legend("topright", legend = c("Média Amostral", "Limite Crítico"), col = c("red", "blue"), lwd = 2, lty = 2)
cat("Comparando a probabilidade empírica de erro tipo I (0.0506) com o nível de significância especificado
(0.05),verificamos que o teste está fornecendo resultados consistentes com o nível de significância desejado.
Se a probabilidade empírica de erro tipo I estiver próxima do nível de significância especificado,
isso sugere que o teste está funcionando conforme esperado, e a taxa de erro tipo I observada
está em linha com a taxa de erro tipo I desejada.")## Comparando a probabilidade empírica de erro tipo I (0.0506) com o nível de significância especificado
## (0.05),verificamos que o teste está fornecendo resultados consistentes com o nível de significância desejado.
## Se a probabilidade empírica de erro tipo I estiver próxima do nível de significância especificado,
## isso sugere que o teste está funcionando conforme esperado, e a taxa de erro tipo I observada
## está em linha com a taxa de erro tipo I desejada.6. Considerando a implementação em R a seguir, descreva detalhadamente sobre o que é, e comente cada linha de instrução do código.
# ----------- CRIANDO DOIS VETORES PARA MAQUINA 1 E MAQUINA 2
maquina1 <- c(15.18, 15.52, 15.42, 14.25, 12.73, 15.51, 10.81, 12.16, 12.85, 12.84,
15.21, 18.20, 11.97, 15.87, 16.80, 13.02, 14.89, 16.65, 14.49, 14.56,
15.62, 14.88, 13.89, 14.72, 18.77, 16.75, 13.17, 12.51, 14.28, 17.66,
15.59, 13.60, 16.76, 14.73, 12.76, 15.92, 18.05, 15.87, 15.38, 13.69,
16.14, 12.86, 11.69, 14.91, 14.93, 19.73, 12.57, 15.34, 16.61, 17.10,
14.98, 13.51, 14.87, 18.88, 15.97, 10.91, 17.85, 16.08, 14.93, 14.96)
maquina2 <- c(16.98, 26.68, 20.38, 18.21, 8.13, 13.73, 9.80, 11.66, 23.74, 10.87,
27.29, 12.14, 24.46, 28.07, 18.21, 13.78, 14.27, 33.99, 7.47, 33.41,
21.60, 12.76, 13.88, 25.31, 9.13, 8.56, 10.31, 20.70, 21.72, 12.35,
25.67, 26.78, 22.39, 14.69, 14.45, 23.34, 16.57, 5.21, 6.67, 11.68,
15.02, 11.48, 20.84, 17.16, 12.46, 18.27, 12.40, 8.54, 6.11, 13.00,
17.13, 15.91, 12.76, 17.60, 21.31, 12.52, 17.69, 20.71, 18.52, 22.92)Criando a função para implementação dos dados das máquinas - comparar as médias de populações (maquina1 x maquinas)
B <- 2000 # ------- Criando uma amostra de 2000 numeros de repetições para o bootstrap
dist1 <- dist2 <- numeric(length = B) # Criando dois vetores VAZIOS de B, também poderia ser feito dist1 <- dist2 NULL
for(i in 1:B){ # Varrendo todo intervalo definido acima em B , ou seja um LOOP repetindo B VEZES
amostra1 <- sample(maquina1, replace = TRUE) # Criando repetições com REPOSIÇÃO de numeros da maquina1
amostra2 <- sample(maquina2, replace = TRUE) # Criando repetições com REPOSIÇÃO de numeros da maquina2
dist1[i] <- mean(amostra2) - mean(amostra1) # Calculando a média das DIFERENÇAS entre amostra2 x amostra1 e guarda no objeto dist1
dist2[i] <- sd(amostra2)/sd(amostra1) # Calculando o DESVIO PADRÃO entre amostra2 x amostra1 e guarda no objeto dist2
}
# Mostrando os 10 primeiros sorteios de cada objeto
head(dist1, 10)## [1] 1.116500 1.504000 0.283000 2.363500 3.520667 2.010167 0.957000 1.039667
## [9] 3.228500 3.478000head(dist2, 10)## [1] 3.440943 3.437329 3.154347 3.712532 3.219498 2.879831 3.444998 3.133958
## [9] 3.354357 3.670715Calculando os quantis das diferenças das médias e quantis do desvio padrão , com nível de 95%
alpha <- 0.05 # tamanho do NÍVEL DE SIGNIFICÂNCIA para os intervalos de confiança
quantile(dist1, c(alpha/2, 1 - alpha/2)) # Calcula os quantis da dist. das diferenças de medias dist1 p/ obter um intervalo de confiança de 95%## 2.5% 97.5%
## 0.109950 3.609525quantile(dist2, c(alpha/2, 1 - alpha/2)) # Calcula os quantis da dist. do desvio padrão de dist2 p/ obter um intervalo de confiança de 95% dos desvios padrão## 2.5% 97.5%
## 2.665072 4.356087# "As médias populcionais são diferentes pq se o zero pertencesse ao intervalo as médias populacionais seriam iguais, podemos perceber que não" 7.O conjunto de dados ‘DAAG::ironslag’ tem 53 observações de medições em pares do teor de ferro obtidos por dois métodos: químico e magnético. A ajuste o modelo linear chemical ∼ magnetic e obtenha os intervalos de confiança bootstrap para os coeficientes do modelo (β0 e β1). As distribuições de bootstrap dos estimadores dos parâmetros β0 e β1 são obtidas seguindo os seguintes passos: Pares (xi, yi), i = 1, . . . , n são os dados amostrais observados. Para r = 1 até R:
data("ironslag")
datatable(ironslag, extensions = "Buttons",
options = list(dom = "Bfrtip",
buttons = c("pdf"), pageLength = 5))modelo linear quimico em relação ao magnético
modelo <- lm(chemical ~ magnetic, data = ironslag)
modelo##
## Call:
## lm(formula = chemical ~ magnetic, data = ironslag)
##
## Coefficients:
## (Intercept) magnetic
## 8.9565 0.5866summary(modelo)##
## Call:
## lm(formula = chemical ~ magnetic, data = ironslag)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.5828 -2.6893 -0.3825 2.7240 6.6572
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8.95650 1.65235 5.420 1.63e-06 ***
## magnetic 0.58664 0.07624 7.695 4.38e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.464 on 51 degrees of freedom
## Multiple R-squared: 0.5372, Adjusted R-squared: 0.5282
## F-statistic: 59.21 on 1 and 51 DF, p-value: 4.375e-10Coeficientes do modelo dos dois coeficientes:
**chemical *intercept variável dependente (x) e magnetic** variável independente (y)
coef(modelo) ## (Intercept) magnetic
## 8.9565011 0.5866413coef_mod <- data.frame(coef(modelo)) # criando um dataframe para facilitar usar o cat
cat("Em média, para cada unidade que 'magnetic' que se aumenta o 'chemical' aumenta cerca de:", coef_mod$coef.modelo.[2], "\n", "unidades, supondo que outras variáveis permaneçam constantes")## Em média, para cada unidade que 'magnetic' que se aumenta o 'chemical' aumenta cerca de: 0.5866413
## unidades, supondo que outras variáveis permaneçam constantes\(y = \beta0 + \beta1x\) modelo de regressão para encontrar o intervalo de confiança de \(\beta0\) e \(\beta1\)
x <- ironslag$magnetic
y <- ironslag$chemical
n <- dim(ironslag)[1]
beta0 <- beta1 <- numeric(1000)
for(i in 1:1000){
indice <- sample(1:n, replace = T) # Sorteie aleatoriamente n inteiros {i∗1, . . . , i∗n} de {1, . . . , n} com reposição
# Defina x∗j = xi∗je y∗j = y∗j, j = 1, . . . , n.
novo_x <- x[indice]
novo_y <- y[indice]
# Ajuste um modelo de regressão de mínimos quadrados ordinários aos pares
# (x∗j, y∗j), j = 1, . . . , n para obterestimativas βˆ(r)∗0e βˆ(r)∗1e o erro padrão residual s∗r para a r-ésima amostra.
novo_mod <- lm(novo_y ~ novo_x)
beta0[i] <- coef(novo_mod)[1]
beta1[i] <- coef(novo_mod)[2]
}
betaZero <- quantile(beta0, c(0.025, 0.975)) # Com 95% de confiança o parametro beta0
cat("Com 95% de confiança o paramentro beta0 está entre",betaZero, "\n")## Com 95% de confiança o paramentro beta0 está entre 5.273592 11.78082# Intervalos de confiança para os coeficientes do modelo de regressão linear
confint(modelo)## 2.5 % 97.5 %
## (Intercept) 5.6392798 12.2737223
## magnetic 0.4335801 0.7397025