Ist das eine Monte-Carlo-Simulation?

Vergleichen wir also zwei Normalverteilungen

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

Die mc.p.-Werte, die unter unser Alpha (0,05) geteilt durch x fallen, würden dann die Fehlerrate vom Typ 1 ergeben. Unser H0 ist a.samples> = b.samples. (Inspiriert von https://www.countbayesie.com/blog/2015/3/3/6-amazing-trick-with-monte-carlo-simulations )

Ich dachte jedoch, dass eine Montecarlo-Simulation die folgenden Schritte ausführen muss:

Algorithmus:

1. Richten Sie eine Verteilung für die Daten f () oder f (θ) und einige H0 ein
2. Wiederholen Sie die folgenden zwei Schritte viele Male: (a) Simulieren Sie einen Datensatz gemäß H0. (B) Berechnen Sie T (x) anhand der simulierten Daten
3. Addiere T (X), das aus den Probendaten ausgewertet wurde
4. Bestellen Sie alle T (x) s
5. Der p-Wert ist der Anteil der T (x) s, der extrem oder extremer ist als der aus den Probendaten

Daher ist das erste Code-Snippet keine echte Monte-Carlo-Simulation? und ist der p-Wert gültig, denn wenn Sie ihn grafisch darstellen, erhalten Sie nicht die erwartete Fehlerrate von 5% Typ 1, die man für einen statistischen Test erwarten könnte.

Dies ist eine Monte-Carlo-Simulation, obwohl ich bezweifle, dass sie das tut, was Sie wollen, und wirklich keinen praktischen Nutzen hat. Sie vergleichen 10000 Einzelprobenstudien und bestimmen, wie viele dieser einzelnen beobachteten Werte im Durchschnitt höher sind. Daher ist es wahrscheinlich besser, Ihren Code als den folgenden weniger effizienten Code zu konzipieren:

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)

Der obige Code sollte dies zeigen, wenn die Verteilungen gleich sind $a$ ist höher als $b$ 50% der Zeit, aber diese Art von Ergebnis ist in der Praxis im Wesentlichen nutzlos, weil Sie nur haben würden $N=2$und statistische Schlussfolgerungen sind nicht anwendbar (dh es gibt keine Varianz innerhalb jeder Gruppe, die zur Quantifizierung der Stichprobenunsicherheit verfügbar ist).

Was Sie vermissen, ist der Vergleich zweier interessierender zusammenfassender Statistiken , um deren Stichprobeneigenschaften zu bestimmen. Dies ist im Allgemeinen das Thema der statistischen Inferenz und erfordert mindestens eine Mindestanzahl von Datenpunkten, um eine Form der Stichprobenunsicherheit zu quantifizieren.

Derzeit ist dies in der Regel ein unabhängiger Standard $t$-Versuchsaufbau. Sie können also eine Reihe von Dingen vergleichen, z. B. wie oft der erste Gruppenmittelwert höher ist als der zweite, wie oft a$t$-Test gibt a $p < \alpha$ Ergebnis (oder analog, wie oft die $|t|$ Verhältnis ist größer als der Grenzwert) und so weiter.

ZB wenn $n=20$ für jede Gruppe und die Verteilungen sind dann in der Bevölkerung gleich

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)

Das Herumspielen mit den anderen Parametern, z. B. das Ändern des ersten Gruppenmittelwerts auf 1, ändert die Art der Simulation (Ändern von einer Fehlersimulation vom Typ I, wie sie jetzt ist, zu einer Leistungsstudie).