normaltesting.R

Jim — Nov 2, 2013, 10:12 AM

# Testing World Series example

# Let mu be the length (minutes) of current World Series games

# Prior beliefs:
# Look at lengths of all games in 2012 season
# Mean length = 180.2 minutes
# If I believe these are similar in length to lengths of WS series,
# my prior is N(180, 15)

# Observe actual lengths of 2013 WS games:
# From baseball-reference.com
# 3:17, 3:05, 3:54, 3:34, 2:52, 3:15
# assume sampling sd is 27

hours <- c(3, 3, 3, 3, 2, 3)
minutes <- c(17, 5, 54, 34, 52, 15)
time <- 60 * hours + minutes

library(LearnBayes)

# illustration of Bayesian one-sided test

m0 <- 180                                # testing mu <= m0
normpar <- c(180, 15)                    # mean and sd of normal prior
data <- c(mean(time), length(time), 27)  # ybar, n, sigma

mnormt.onesided(m0, normpar, data)
$BF
[1] 0.08342

$prior.odds
[1] 1

$post.odds
[1] 0.08342

$postH
[1] 0.077

# illustration of Bayesian two-sided test

m0 <- 180                                # testing mu = m0
prior.prob <- 0.5                        # prior prob of hypothesis
tau <- 10                                # sd of normal prior under alt
data <- c(mean(time), length(time), 27)  # ybar, n, sigma

mnormt.twosided(m0, prior.prob, tau, data)
$bf
[1] 0.6662

$post
[1] 0.3998

# compare with frequentist p-value

z <- (mean(time) - 180) * sqrt(length(time)) / 27
p.value <- 2 * (1 - pnorm(z))
p.value
[1] 0.07688

# check the sensitivity of posterior prob of H with respect to
# prior standard deviation

m0 <- 180                                # testing mu = m0
prior.prob <- 0.5                        # prior prob of hypothesis
data <- c(mean(time), length(time), 27)  # ybar, n, sigma

myfunction <- function(tau)
  mnormt.twosided(m0, prior.prob, tau, data)$post

curve(myfunction(x), 0.1, 40,
      xlab = "Prior Standard Deviation",
      ylab = "Posterior Prob of H")

plot of chunk unnamed-chunk-1


# we see that P(H | data ) >= .38 for all tau