Este codigo tem por objetivo realizar simulações pelo metodo de Monte Carlo com intuito de verificar o poder do teste de algumas distribuições.

Teste t

Para todos os casos utilizaremos um nivel de significância de \(\alpha= 0,05\) Como tamanho das amostras foi definido: 10,20,100 e 200. E para cada distribuição um valor experado de \(\mu\) foi definido.

Monte Carlo

Para a simulação de Monte Carlos seram feitas 1000 repetições.

Distribuição \(U_{[0,1]}\)

Para a distribuição uniforme utilizaremos \(\mu = \frac{1}{2}\)

alpha <- 0.05
n <- c(10,20,100,200)
MC <- 1000 

mu0 <- seq(0.1, 0.9, by = 0.05)

resultado <- array(NA, dim = c(MC, length(mu0), length(n)))

for(k in 1:length(mu0)){
  for(j in 1:length(n)){
    for(i in 1:MC){
      amostra <- runif(n[j],min = 0, max = 1)
      resultado[i,k,j] <- t.test(amostra, mu = mu0[k],
                                 alternative = "less")$p.value < alpha
    }
  }
}

N10 <- cbind(mu0, colMeans(resultado[,,1]))
N20 <- cbind(mu0, colMeans(resultado[,,2]))
N100 <- cbind(mu0, colMeans(resultado[,,3]))
N200 <- cbind(mu0, colMeans(resultado[,,4]))

plot(N10, type = "l", axes = F, lwd = 2, col= 'chocolate3',
     xlab = "Valores de mu0", ylab = "Poder do teste", main =  "Distribuição Uniforme")
lines(N20, type = "l", lwd = 2, col= 'blue')
lines(N100, type = "l", lwd = 2, col= 'green')
lines(N200, type = "l", lwd = 2, col= 'gold2')
abline(v = 0.5, lwd = 3, col= 'red')
abline(h = 0.05, lwd = 2, col = "grey80", lty = 2)
axis(1, seq(0.1, 0.9, .1), lwd = 2)
axis(2, seq(0.1, 0.9, .2), lwd = 2)
box(lwd = 3)

legend("topleft", ncol = 1, col = c('chocolate3', 'blue', 'green', 'gold2'),
       legend = c('N = 10', 'N = 20', 'N = 100', 'N = 200'),
       bty="n", lty = 1, lwd = 2)

Distribui??o \(Exp(1/2)\)

Para a distribuição exponencial utilizaremos \(\mu = \frac{1}{2}\).

alpha <- 0.05         
n <- c(10,20,100,200) 
MC <- 1000            
mu0 <- seq(-1, 1, by = 0.025)

resultado <- array(NA, dim = c(MC, length(mu0), length(n)))

for(k in 1:length(mu0)){
  for(j in 1:length(n)){
    for(i in 1:MC){
      amostra <- rexp(n[j], 2)
      resultado[i,k,j] <- t.test(amostra, mu = mu0[k],
                                 alternative = "less")$p.value < alpha
    }
  }
}

N10 <- cbind(mu0, colMeans(resultado[,,1]))
N20 <- cbind(mu0, colMeans(resultado[,,2]))
N100 <- cbind(mu0, colMeans(resultado[,,3]))
N200 <- cbind(mu0, colMeans(resultado[,,4]))

plot(N10, type = "l", axes = F, lwd = 2, col= 'chocolate3',
     xlab = "Valores de mu0", ylab = "Poder do teste",
     xlim = c(-0.1,1), ylim = c(0,1),main =  "Distribuição Exponencial")
lines(N20, type = "l", lwd = 2, col= 'blue')
lines(N100, type = "l", lwd = 2, col= 'green')
lines(N200, type = "l", lwd = 2, col= 'gold2')
abline(v = 0.5, lwd = 3, col= 'red')
abline(h = 0.05, lwd = 2, col = "grey80", lty = 2)
axis(1, seq(-.1, 1, .2), lwd = 2)
axis(2, seq(0.1, 0.9, .2), lwd = 2)
box(lwd = 3)

legend("topleft", ncol = 1, col = c('chocolate3', 'blue', 'green', 'gold2'),
       legend = c('N = 10', 'N = 20', 'N = 100', 'N = 200'),
       bty="n", lty = 1, lwd = 2)

Distribuição \(N(1/2,1)\)

Para a distribuição normal utilizaremos \(\mu = \frac{1}{2}\) e \(\sigma= 1\), que são estimados por \(\overline{X}\).

alpha <- 0.05         
n <- c(10,20,100,200) 
MC <- 1000            
mu0 <- seq(-1, 3, by = 0.025)

resultado <- array(NA, dim = c(MC, length(mu0), length(n)))

for(k in 1:length(mu0)){
  for(j in 1:length(n)){
    for(i in 1:MC){
      amostra <- rnorm(n[j], mean = 0.5)
      resultado[i,k,j] <- t.test(amostra, mu = mu0[k],
                                 alternative = "less")$p.value < alpha
    }
  }
}

N10 <- cbind(mu0, colMeans(resultado[,,1]))
N20 <- cbind(mu0, colMeans(resultado[,,2]))
N100 <- cbind(mu0, colMeans(resultado[,,3]))
N200 <- cbind(mu0, colMeans(resultado[,,4]))

plot(N10, type = "l", axes = F, lwd = 2, col= 'chocolate3',
     xlab = "Valores de mu0", ylab = "Poder do teste",
     xlim = c(-1,1), ylim = c(0,1), main =  "Distribuição Normal")
lines(N20, type = "l", lwd = 2, col= 'blue')
lines(N100, type = "l", lwd = 2, col= 'green')
lines(N200, type = "l", lwd = 2, col= 'gold2')
abline(v = 0.5, lwd = 3, col= 'red')
abline(h = 0.05, lwd = 2, col = "grey80", lty = 2)
axis(1, seq(-1.1, 1, .2), lwd = 2)
axis(2, seq(0.1, 0.9, .2), lwd = 2)
box(lwd = 3)

legend("topleft", ncol = 1, col = c('chocolate3', 'blue', 'green', 'gold2'),
       legend = c('N = 10', 'N = 20', 'N = 100', 'N = 200'),
       bty="n", lty = 1, lwd = 2)

Distribuição \(Cauchy(1/2,1)\)

Para a distribuição de Cauchy utilizaremos como parametros \(\alpha=1/2\) e \(\beta=1\), pois a distribuição nao apresenta \(\mu\).

alpha <- 0.05
n <- c(10,20,100,200) 
MC <- 1000
mu0 <- seq(-3.5, 5.5, by = 0.5)

resultado <- array(NA, dim = c(MC, length(mu0), length(n)))

for(k in 1:length(mu0)){
  for(j in 1:length(n)){
    for(i in 1:MC){
      amostra <- rcauchy(n[j], location = 0.5, scale = 1)
      resultado[i,k,j] <- t.test(amostra, mu = mu0[k],
                                 alternative = "less")$p.value < alpha
    }
  }
}

N10 <- cbind(mu0, colMeans(resultado[,,1]))
N20 <- cbind(mu0, colMeans(resultado[,,2]))
N100 <- cbind(mu0, colMeans(resultado[,,3]))
N200 <- cbind(mu0, colMeans(resultado[,,4]))

plot(N10, type = "l", axes = F, lwd = 2, col= 'chocolate3',
     xlab = "Valores de mu0", ylab = "Poder do teste", main =  "Distribuição Cauchy")
lines(N20, type = "l", lwd = 2, col= 'blue')
lines(N100, type = "l", lwd = 2, col= 'green')
lines(N200, type = "l", lwd = 2, col= 'gold2')
abline(v = 0.5, lwd = 3, col= 'red')
abline(h = 0.05, lwd = 2, col = "grey80", lty = 2)
axis(1, seq(-8, 8, .5), lwd = 2)
axis(2, seq(0.1, 0.9, .2), lwd = 2)
box(lwd = 3)

legend("topleft", ncol = 1, col = c('chocolate3', 'blue', 'green', 'gold2'),
       legend = c('N = 10', 'N = 20', 'N = 100', 'N = 200'),
       bty="n", lty = 1, lwd = 2)

Distribuição \(LogN(1/2,1)\)

Para a distribuição normal utilizaremos \(\mu = \frac{1}{2}\) e \(\sigma= 1\), que são estimados por \(\overline{X}\).

alpha <- 0.05         
n <- c(10,20,100,200) 
MC <- 1000           
mu0 <- seq(0, 3, by = 0.05)
resultado <- array(NA, dim = c(MC, length(mu0), length(n)))
for(k in 1:length(mu0)){
  for(j in 1:length(n)){
    for(i in 1:MC){
      amostra <- rlnorm(n[j], meanlog = 0.5, sdlog = 1)
      resultado[i,k,j] <- t.test(amostra, mu = mu0[k],
                                 alternative = "less")$p.value < alpha
    }
  }
}
N10 <- cbind(mu0, colMeans(resultado[,,1]))
N20 <- cbind(mu0, colMeans(resultado[,,2]))
N100 <- cbind(mu0, colMeans(resultado[,,3]))
N200 <- cbind(mu0, colMeans(resultado[,,4]))
plot(N10, type = "l", axes = F, lwd = 2, col= 'chocolate3',
     xlab = "Valores de mu0", ylab = "Poder do teste",
     main =  "Distribuição Log Normal")


lines(N20, type = "l", lwd = 2, col= 'blue')
lines(N100, type = "l", lwd = 2, col= 'green')
lines(N200, type = "l", lwd = 2, col= 'gold2')
abline(v = 0.5, lwd = 3, col= 'red')
abline(h = 0.05, lwd = 2, col = "grey80", lty = 2)
axis(1, seq(-1, 3, .5), lwd = 2)
axis(2, seq(0, 0.9, .05), lwd = 2)
box(lwd = 3)
legend("topleft", ncol = 1, col = c('chocolate3', 'blue', 'green', 'gold2'),
       legend = c('N = 10', 'N = 20', 'N = 100', 'N = 200'),
       bty="n", lty = 1, lwd = 2)