Illustrates Model-Data Simulation to Learn About a Player's Batting Ability Based on a "ofer" Slump

model.data.sim.R

# Script to Learn About Ryan Howard's Batting Ability from a "0 for 35" Slump

# Uses a function from the BayesTestStreak package
# install_github("bayesball/BayesTestStreak")

library(MASS)
library(BayesTestStreak)


# Simulate 500 at-bats with a constant hitting probability p = 0.250 
# -- find the longest ofer.

one.sim <- function(p=.25, n=500){
require(BayesTestStreak)
max(find.spacings(rbinom(n, size=1, prob=p))$y)
}

# Repeat this simulation 1000 iterations 
# graph the maximum ofers and overlay the observed value.

S <- replicate(1000, one.sim())
truehist(S,
main="Longest Ofer in 500 AB if p = 0.250")
abline(v=35, lwd=3, col="red")
text(30, .07, "OBSERVED", col="red")

# Here is a different question.  
# What can we learn about Howard's true batting average on the 
# basis of the data "oh for 35 or greater"?

# Illustrate model/data simulation.

# Suppose Howard's true batting p is any value between 
# p=.10 and p=.40 and we assign each value the same prior probability.

# Here is the prior:

p <- seq(.1, .4, .01)
prior <- 1/length(p) + 0 * p
plot(p, prior, type="h", col="red", 
lwd=6, main="The Prior on Howard's True AVG")

model.data.sim <- function(){
require(BayesTestStreak)
p <- sample(seq(.1, .4, .01), size=1)
x <- max(find.spacings(rbinom(500, size=1, prob=p))$y)
c(p, x)
}

many.sim <- function(ofer, iter=100000){
S <- replicate(iter, model.data.sim())
S.p <- S[1, S[2, ] >= ofer]
truehist(S.p, 
main=paste("Howard's Posterior Given Ofer =", ofer))
}


# What is a player goes "0 for 15+" in a 500-AB season?

many.sim(15)

# What if a player goes "0 for 35+" in a 500-AB season?

many.sim(35)