Binomial testing with buttered toast
Rasmus' post of last week on binomial testing made me think about p-values and testing again. In my head I was tossing coins, thinking about gender diversity and toast. The toast and tossing a buttered toast in particular was the most helpful thought experiment, as I didn't have a fixed opinion on the probabilities for a toast to land on either side. I have yet to carry out some real experiments.Suppose I tossed 6 buttered toasts and observed that all but one toast landed on their buttered side.
Now I have two questions:
- Would it be reasonable to assume that there is a 50/50 chance for a toast to land on either side?
- Which probability should I assume?
The probability of observing one ore more B (right tail event) is:
(gt <- 1-1/2^6)
# [1] 0.984375
and the probability of observing one or fewer B (left tail event) is:(lt <- 1/2^6*choose(6,1) + 1/2^6)
# [1] 0.109375
while the probability of either extreme event, one or fewer B (or U), is:2*min(c(gt, lt))
# [1] 0.21875
In summary, if the toast has an equally probability to fall on either side, then there is 22% chance to observe one or fewer B (or U) in 6 tosses. That's not that unlikely and hence I would not dismiss the hypothesis that the toast is fair, despite the fact that the sample frequency is only 1/6.Actually, the probabilities I calculated above are exactly the p-values I get from the classical binomal tests:
## Right tail event
binom.test(1, 6, alternative="greater")
## Left tail event
binom.test(1, 6, alternative="less")
## Double tail event
binom.test(1, 6, alternative="two.sided")
Additionally I can read from the tests that my assumption of a 50% probability is on the higher end of the 95 percent confidence interval. Thus, wouldn't it make sense to update my belief about the toast following my observations? In particular, I am not convinced that a 50/50 probability is a good assumption to start with. Arguably the toast is biased by the butter.Here the concept of a conjugate prior becomes handy again. The idea is to assume that the parameter \(\theta\) of the binomial distribution is a random variable itself. Suppose I have no prior knowledge about the true probability of the toast falling on either side, then a flat prior, such as the Uniform distribution would be reasonable. However, the beta distribution with parameter \(\alpha=1\) and \(\beta=1\) has the same property and is a conjugate to the binomial distribution with parameter \(\theta\). That means there is an analytical solution, in this case the posterior distribution is beta-binomial with hyperparaemters:
\[\alpha':=\alpha + \sum_{i=1}^n x_i,\; \beta':=\beta + n - \sum_{i=1}^n x_i,\]
and the posterior predictor for one trial is given as
\[\frac{\alpha'}{\alpha' + \beta'}\]
so in my case:
alpha <- 1; beta <- 1; n <- 6; success <- 1
alpha1 <- alpha + success
beta1 <- beta + n - success
(theta <- alpha1 / ( alpha1 + beta1))
# [1] 0.25
My updated believe about the toast landing on the unbuttered side is a probability of 25%. That's lower than my prior of 50% but still higher than the sample frequency of 1/6. If I would have more toasts I could run more experiments and update my posterior predictor.I get the same answer from Rasmus'
bayes.binom.test
function:> library(BayesianFirstAid)
> bayes.binom.test(1, 6)
Bayesian first aid binomial test
data: 1 and 6
number of successes = 1, number of trials = 6
Estimated relative frequency of success:
0.25
95% credible interval:
0.014 0.527
The relative frequency of success is more than 0.5 by a probability of 0.061
and less than 0.5 by a probability of 0.939
Of course I could change my view on the prior and come to a different conclusion. I could follow the Wikipedia article on buttered toast and believe that the chance of the toast landing on the buttered side is 62%. I further have to express my uncertainty, say a standard deviation of 10%, that is a variance of 1%. With that information I can update my belief of the toast landing on the unbuttered side following my observations (and transforming the variables):x <- 0.38
v <- 0.01
alpha <- x*(x*(1-x)/v-1)
beta <- (1-x)*(x*(1-x)/v-1)
alpha1 <- alpha + success
beta1 <- beta + n - success
(theta <- alpha1 / ( alpha1 + beta1))
# [1] 0.3351821
I would conclude that for my toasts / tossing technique the portability is 34% to land on the unbuttered side. In summary, although their is no sufficient evidence to reject the hypothesis that the 'true' probability is not 50% (at the typical 5% level), I would work with 34% until I have more data. Toast and butter please!
Session Info
R version 3.0.1 (2013-05-16)
Platform: x86_64-apple-darwin10.8.0 (64-bit)
locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods
[7] base
other attached packages:
[1] BayesianFirstAid_0.1 rjags_3-12 coda_0.16-1
[4] lattice_0.20-24
loaded via a namespace (and not attached):
[1] grid_3.0.1 tools_3.0.1
7 comments :
I've also been thinking of good example of binary "processes" which are not 50%/50%, coins are such a bad example. Toast is good example!
I (almost) found some data on this. Though it does not mentions the exact numbers it seems like out of almost 10000 trials around 6200 toasts fell with the buttered side down: http://www.telegraph.co.uk/news/uknews/1331810/Breakfast-at-Murphys-or-why-the-toast-lands-butter-side-down.html. For quick beta stuff I often use the rbeta function in R that samples from a beta distribution. This would, for example, result in the median posterior and a 95% equitail credible interval:
quantile( rbeta(99999, 1 + 1, 1 + 6), c(0.05, 0.5, 0.975)
Thanks Rasmus. That's the article the Wikipedia is referring to.
But why simulate when you can have the answers directly via the qbeta function:
qbeta(c(0.025, 0.5, 0.975), 1 + 1, 1 + 6)
True, true. I guess I'm more of a random guy. Livin' on the edge, YOLO, and that kind of stuff.
regarding illustrative binary processes that aren't 50:50, what about the thumbtack flipping experiment, Pr(land on flat side). I saw that somewhere ....
Hello Markus,
I'm trying to find the optimum values (a and b) for a function's minumum value or a constant minimum level. I used the optim and optimize function but I can't find the optimal pairs of the parameters of a and b. My function is some more complicated and I create this function as min.RSS(a)(b) (two functions are nested) when I used the optim function and I received the following error from this code: error in optim(c(1,1), min.RSS) : cannot coerce type 'closure' to vector of type 'double'
Could you please help me on finding the optimal pairs of parameters.
Thanks in advance
Hi Jenny,
You will have to provide a minimal reproducible example for anyone to help you. Those questions are best ask by email, or at stack overflow or R-help.
Hi Sir,
Could you please tell me how to use optim for finding the maximum of a function?
Thanks
Anand
Post a Comment