###### Return and plot the posterior distribution of the number of
###### dice rolled for a given sample of rolls and prior beliefs.
##
#### Arguments
##
## rolls = vector of roll totals, assuming they are independent and
## made with the same number of dice
##
## prior.min = smallest plausible number of dice (prior to seeing the
## actual rolls)
##
## prior.max = greatest plausible number of dice (prior to seeing the
## actual rolls)
##
## prior.p = a vector of prior probabilities assigned to each value
## between prior.min and prior.max (inclusive); defaults to uniform
##
## plot = plot the posterior mass?
##
## cutoff = don't bother plotting possibilities with posterior
## probability less than this
##
#### Caveats
##
## * Calculating probabilities for sums of dice is a hassle, so I am
## using a Normal approximation here. It should work fine if
## prior.min > 3 or so, and definitely if prior.min > 6.
##
## * Rolls are assumed to be,
##
## 1) Independent: Probably true for dice rolls.
##
## 2) Identically distributed: The true underlying number of dice is
## the same for each roll, so no hidden challenge dice or boosts.
##
## 3) Missing at random for a given d: That is, rolls that came out
## too low or too high are not hidden from us, though it's OK if
## the QMs didn't bother to roll because the fight was too uneven
## or probabilities of all but one outcome was negligible.
##
d.post <- function(rolls, prior.min, prior.max, prior.p=1,
plot=TRUE, cutoff=.005){
prior.p <- rep(prior.p, length.out=prior.max-prior.min+1)
prior <- function(d) prior.p[d-prior.min+1]
## Normal approximation: should work well enough for d>3. I probably
## could have just used Normal PDF here, but this should be a tiny
## bit more accurate.
lik <- function(d){
m <- 50.5*d
s <- sqrt(d*(100^2-1)/12)
prod((pnorm(rolls+.5, m, s)
-pnorm(rolls-.5, m, s)))/
(pnorm(d*100+.5, m, s)
-pnorm(d-.5, m, s))^length(rolls)
}
post <- sapply(prior.min:prior.max, function(d) prior(d)*lik(d))
names(post) <- prior.min:prior.max
post <- post/sum(post)
if(plot)
plot((prior.min:prior.max)[post>=cutoff],
post[post>=cutoff],
xlab="Dice", ylab="Posterior probability",
type="h")
post
}