Let's say I want to make a football simulator based on real-life data.
Say I have a player who averages 5.3 yards per carry with a SD of 1.7 yards.
I'd like to generate a random variable that simulates the next few plays. eg: 5.7, 4.9, 5.3, etc.
What stats terms to I need to look up to pursue this idea? Density function? The normal curve estimates what boundaries the data generally fall within, but how do I translate that into simulation of subsequent data points?
Thanks for any guidance!
Of course you can use rnorm() in R, but it may be easier to understand how drawing from a pdf works by using the probability integral transform.
Basically, once we specify the structure of the pdf, we can transform this into a cdf (empirically, to ignore what the equation is), and because the values of the cdf have unique values from 0 to 1, we can back-calculate a draw from the original pdf by matching random draws from 0 to 1, with the cdf.
This way, you only need to have a RNG from 0 to 1, and the function of the pdf, and you're set. Here is the R code:
x <- seq(-4, 4, len = 1000)
f <- function(x, mu = 0, sigma = 1) {
out <- 1 / sqrt(2*pi*sigma^2) * exp(-(x - mu)^2 / (2*sigma^2))
out
}
x.ecdf <- cumsum(f(x)) / sum(f(x))
out <- vector()
y <- runif(100)
for (i in 1:length(y)) {
out[i] <- which((y[i] - x.ecdf)^2 == min((y[i] - x.ecdf)^2))
}
par(mfrow = c(1,2))
plot(x, x.ecdf)
hist(x[out], breaks = 20)
alt text http://probabilitynotes.files.wordpress.com/2010/08/rnormish.png
A player's yardage is unlikely to be anywhere near normally distributed. If it were, your guy averaging 5.3 give or take 1.7 yards would almost never lose yards or gain more than 11 yards on any play in the entire season. Gone is the excitement of the game, to be replaced by some statistical mediocrity. If football were played like this, a team's chances of making a set of downs would be almost certain; there would almost never be a loss of downs; and the game would simply be determined by who won the initial coin flip and got on the field first.
Why not just draw a value at random from a list of the player's recent gains (and losses)? It's fairly easy to program: you just have to generate a uniformly distributed integer to index into an array of the gains. It doesn't require any kind of statistical model--no need to fit anything. It can account for change in the player's ability over time (just by selecting which time period you will use to draw the data from). And it's obviously driven by "real-live data."
You need a random number generator for the standard normal. Either you can supply mean and standard deviation as arguments to the function, or you simply scale yourself by multiplying with the latter for the variability and adding the former for for the central location.
Here is a quick example of the former approach:
> set.seed(42)
> x <- rnorm(1000, 5.3, 1.7) # 1000 draws of N(5.5, 1.7)
> print(c(mean=mean(x), sd=sd(x)))
mean sd
5.2561 1.7043
>
If you want a realistic simulation you need to find a distribution that describes the real process good enough (a model).
When a real player makes a move he will on average (e.g.) throw X yards, with a standard deviation of Y. This does however not mean that the distribution of throws is a normal distribution. You should plot the (histogrammed) throw distribution and plot your fit (a Normal distribution with mean X and σ=Y) and determine if it fits good enough. If not find a distribution that describes the real data better.
Once you have that down you need to generate random numbers from the distribution you determined.
EDIT
If you data is complete enough you could maybe avoid having to create a complete model. You would create a frequency distribution of your data (a histogram) and then do rejection sampling directly from it.
if you decide to generate your distribution from the data you have observed, your model will never spit out a "tail value", ie, something outside the range of what you have observed.
your example data: average 5.3 yards per carry with a SD of 1.7 yards
will have a max and a min, say 10 and 2. in that case, your calculated distribution will not have any weight in the tails, and will be unable to generate a value of 11 or 1. maybe this is ok, but it prevents you from ever generating one of those super-human events that everyone loves to watch.
the functional form for the standard normal has tails that go to infinity, so if you assume the distribution is normal, you will be able to generate simulated values (at very small probabilities) that are higher or lower than your observed data.
If you have all the relevant data rather than just summary data such as mean, SD etc. you could create your own distribution model from the real life data you have. Sort the data (y values) from lowest to highest and equally space them between 0 and 1 (x values). Then solve to find the coefficients of an nth order polynomial curve fit to the data (or several in piecewise fashion over parts of the data if necessary). Once you have these, it would be a simple matter to use a uniform random number generator to generate an x value between 0 and 1, and to plug this value into the polynomial equation to get a random y value that would approximate a draw from your distribution.
