Esta é uma simulação de Monte Carlo?

Então, vamos comparar duas distribuições normais

Do this x times: 

runs <- 100000
a.samples <- rnorm(runs, mean = 5) 
b.samples <- rbeta(runs, mean = 0) 
mc.p.value <- sum(a.samples > b.samples)/runs

Os valores mc.p. abaixo de nosso alfa (0,05) dividido por x dariam a taxa de erro do tipo 1. Nosso H0 é a.samples> = b.samples. (Inspirado em https://www.countbayesie.com/blog/2015/3/3/6-amazing-trick-with-monte-carlo-simulations )

Mas, pensei que uma simulação montecarlo tivesse que seguir os seguintes passos:

Algoritmo:

1. Configure alguma distribuição para os dados, f () ou f (θ) e alguns H0
2. Repita as duas etapas a seguir várias vezes: (a) Simule um conjunto de dados de acordo com H0 (b) Calcule T (x) usando os dados simulados
3. Adicione T (X) avaliado a partir dos dados de amostra
4. Encomende todos os T (x) s
5. O valor p é a proporção de T (x) s tão extrema ou mais extrema que a dos dados da amostra

Portanto, o primeiro trecho de código não é uma simulação de boa-fé monte carlo? e o valor p é válido porque, se você fizer um gráfico, não obtém a taxa de erro de tipo 1 de 5% esperada que se poderia esperar de um teste estatístico.

Esta é uma simulação de Monte Carlo, embora eu duvide que esteja fazendo o que você deseja e realmente não seja útil. O que você está comparando são 10000 estudos de amostra única e a determinação de quantos desses valores individuais observados são, em média, mais altos. Portanto, provavelmente é melhor conceituar seu código como o seguinte código menos eficiente:

runs <- 10000
result <- logical(runs)

for(i in 1:runs){
  a <- rnorm(1, mean = 0) 
  b <- rnorm(1, mean = 0)
  result[i] <- a > b
} 
mc.p.value <- mean(result)

O código acima, quando as distribuições são iguais, deve mostrar que é maior que 50% do tempo, mas esse tipo de resultado é essencialmente inútil na prática porque você teria apenas e inferências estatísticas não são aplicáveis ​​( isto é, não há variação dentro de cada grupo disponível para quantificar a incerteza da amostragem).$a$$b$$N=2$

O que está faltando é comparar duas estatísticas resumidas de interesse para determinar suas propriedades de amostragem, que geralmente é o tópico de inferência estatística e requer pelo menos um número mínimo de pontos de dados para quantificar alguma forma de incerteza de amostragem.

Tal como está, geralmente é uma configuração independente de teste padrão . Portanto, você pode comparar várias coisas, como quantas vezes a média do primeiro grupo é maior que a segunda, quantas vezes um teste fornece um resultado (ou analogamente, com que frequência a relação é maior do que o ponto de corte) e assim por diante.$t$$t$$p < \alpha$$|t|$

Por exemplo, se para cada grupo e as distribuições forem iguais na população, então$n=20$

runs <- 10000
n <- 20
alpha <- .05
result.p <- result.mean <- numeric(runs)

for(i in 1:runs){
  a <- rnorm(n, mean = 0) 
  b <- rnorm(n, mean = 0)
  result.mean[i] <- mean(a) > mean(b)
  result.p[i] <- t.test(a, b)$p.value
} 
mc.p.value <- mean(result.p < alpha)
mc.a_gt_b.value <- mean(result.mean)

Brincar com outros parâmetros, como alterar a média do primeiro grupo para 1, mudará a natureza da simulação (alterando-a de uma simulação de erro Tipo I, como é agora, para um estudo de potência).