Chapter 5
One way to look at the Monty Hall simulation is to compare the cumulative wins (or cars if you prefer) over a given horizon. In the plot below, I simulated the Monty Hall problem out to 100 trials.
Convinced?
Here is the code:
# clojure
(def stay (five/simulate-monty-hall 1000 true))
(def switch (five/simulate-monty-hall 1000 false))
(def data (map vector (map #(get % :trial) stay) (map #(get % :wins) stay) (map #(get % :wins) switch)))
(def header ["trial" "stay" "switch"])
(util/write-to-csv "tmp/monty-hall-sim.csv" (conj data header))
# R
df = read.csv("tmp/monty-hall-sim.csv", header=TRUE)
ggplot(aes(x=trial, color=variable), data=df) +
geom_line(aes(y=stay,color="stay")) +
geom_line(aes(y=switch,color="switch")) +
scale_colour_manual(values=c("blue", "red")) +
labs(x="Trial", y="Wins", title="Monty Hall Stay vs Switch")
ggsave("plots/monty-hall.png")
I found that with a value of n = 4 the mean loaf weight converged to 1000 as the number of samples increased. For example, after 365 samples, the mean was 1000.82 and the standard deviation was 33.9.
(def r (five/compare-baker-and-poincare 950 50 4 365))
...
:mean 1000.82465,
:stddev 33.90772240333088,
By plotting the CDF of the simulated Poincaré and a normal distribution with the same mean and variance we may be able to determine if shapes of the distributions are different enough to warrant consulting the bread police. Here is the plot:
They look very similar though there are areas above and below the mean where Poincaré distribution is heavier. Based on this alone, I would not notify the bread police. I would want more evidence.
To gather more evidence, we might also perform a t-test on the data to determine if the mean of the baker's distribution is less than that of Poincaré's. Our null hypothesis would be that there is no significant difference between the baker's distribution and Poincaré's. The alternative hypothesis would be that the baker is giving Poincaré heavier loaves.
# In R
t.test(df$bakery, df$poincare, alternative="less")
t = -1.1371, df = 726.093, p-value = 0.1279
alternative hypothesis: true difference in means is less than 0
The p-value returned is not below 0.05. Using an alpha level of 0.94 would lead us to accept the null hypothesis that there is no significant difference between the simulated bakery loafs and those given to Poincaré.
TODO: add the boxplot?
If you go to a dance where partners are paired up randomly, what percentage of opposite sex couples will you see where the woman is taller than the man?
One way to answer this question is to simulate the pairing by sampling N times each from two random variables M and W pairing them together and computing the ratio of count(W > M) / N. As N increases, this ratio should converge to the expected probability that in a random pairing the woman will be taller than the man.
Here is one way to do that:
(let [horizon 1000
sims 100
men-mean 178
men-var 59.4
men-sigma (Math/sqrt men-var)
women-mean 163
woman-var 52.8
women-sigma (Math/sqrt woman-var)
sim (repeatedly sims (fn []
(/ (count (filter
#(> (second %) (first %))
(map vector
(random/sample horizon #(random/normalvariate men-mean men-sigma))
(random/sample horizon #(random/normalvariate women-mean women-sigma)))))
horizon)))]
(float (stats/mean sim)))The resulting value (float (stats/mean sim)) should give us the simulated probability that a randomly selected woman is taller than a randomly selected man.
The value computed was 0.079 or 7.9%.
As I was considering how to simulate this problem I was curious to know if there was an analytical approach to deriving the same probability and confirm the simulated results. It appears that there is. As I was searching I came upon this question on stats.stackexchange.com. The "stress-strength" analysis can be used to determine the probability that a value from one distribution is greater than the value from another. Stated differently, it will compute P(A < B) where A and B are random variables.
The result of the stress-strength equation is a z score. From the z-score we can then compute the p value or the area under the standard normal distribution. Here is the code:
(stats/z->p-value (p/stress-strength-prob men-mean men-var women-mean woman-var))The computed value with the let bindings above is 0.078 or 7.8% which is very close to our simulated value of 7.9%.
From the above, I would conclude that in a little less than 8% of the pairings, the woman will be taller than the man.
At the simplest level this question reduces to the probability of making 10 shots in a row out of 15 when the probability of making each shot is 50%.
This can be simulated with 10^6 trials using (five/streak-sim 15 10 0.5 1000000) which gives a result of ~0.0034.
To confirm the result, I tried determine the likelihood of a 10-in-a-row streak out of 15 tries. The total event space with 15 tries and two possible outcomes is 2^15. To look at the probability of shooting 10 baskets in a row, we need to look at the different ways in which a shooter can make 10 in a row. When we are done, we can divide the events we are interested in by the total possible outcomes to arrive at a probability of those outcomes.
There are several ways in which a shooter could make 10 in a row. For the illustration below I'll us B for a basket and M for a miss and X for either. X will serve as a placeholder in the positions where either a B or a M will satisfy an arrangement with 10 in a row.
BBBBBBBBBBXXXXX => 2^5 ways we can obtain this 10 in a row outcome
MBBBBBBBBBBXXXX => 2^4
XMBBBBBBBBBBXXX => 2 x 2^3
XXMBBBBBBBBBBXX => 2^2 x 2^2
XXXMBBBBBBBBBBX => 2^3 x 2
XXXXXBBBBBBBBBB => 2^4
Adding these up gives us:
2^5 + 2^4 + 2 x 2^3 + 2^2 x 2^2 + 2^3 x 2 + 2^4 = 112, 112/2^15 = 0.00341 which confirms our simulated result ((five/streak-sim 15 10 0.5 1000000)).
We can simulate 1M cohorts and print some of the summary statistics with the following:
(def s (five/cohort-sim 1000000 100 10 1/1000))
(clojure.pprint/pprint (stats/summary s))
{:95th 3,
:stddev 0.9998101239256801,
:mean 0.999154,
:median 1,
:count 1000000,
:max 9,
:75th 2,
:min 0,
:25th 0}
(def cdf (d/cdf s))
(clojure.pprint/pprint cdf)
{0 0.368183,
1 0.736162,
2 0.91991,
3 0.980945, ; <- less than 0.5%
4 0.996322,
5 0.999425,
6 0.999909,
7 0.999991,
8 0.999999,
9 1.0}With the CDF, we are now prepared to answer part two of this question. In a cohort of 100 people over 10 years, how many cases are at or above the 95th percentile? This is another way of asking how many cases represent a p-value (probability by chance alone) of less than 5%.
From the CDF we see that there is no observation at 0.95. We have 2 observations at 0.91 and 3 at 0.98. The number of observations that satisfies a probability of less than 5% is 3 at 0.98 or 2%.
The distributions/ccdf function can also help us determine this value:
(def cdf (d/cdff s :to-float true))
(cdf 0.95 :value)
3For part three, we want to divide the population of 10000 people into 100 cohorts and follow them for 10 years and then answer the question what is the probability that one of the cohorts will have a "statistically significant" cluster? We can use the probability of 3 observations as satisfying the criteria of "statistically significant" with a p-value of 0.019 or 1.9%. To determine the probability that one cluster out of one hundred will meet this criteria we want to determine the probability that one or more clusters will have 3 observations.
One way to do this would be to use the binomial distribution and take the sum of the probabilities of 1, 2, 3...100 clusters having 3 observations.
(reduce + (map #(p/binomial-pmf 100 % (- 1 (cdf 3))) (range 1 101)))
0.853962162354355That's quite surprising. This indicates that there is an 85% probability that 1 or more clusters will have 3 observations.
If our p-value is 1%, we can use the probability of 4 observations.
(reduce + (map #(p/binomial-pmf 100 % (- 1 (cdf 4))) (range 1 101)))
0.30821525755720985This result is also surprising. If we follow 100 cohorts for 10 years, the probability of a "statistically significant" number of observations with a p-value of 1% is 30%.
Another way to do this would be to simulate 100 cohorts of 100 people over 10 years several times, count the number of cohorts with 3 or more (or 4 or more) and average those counts over the simulations.
(float (stats/mean (repeatedly 10000 #(count (filter (fn [x] (>= x 3)) (five/cohort-sim 100 100 10 1/1000))))))
8.012
(float (stats/mean (repeatedly 1000 #(count (filter (fn [x] (>= x 4)) (five/cohort-sim 100 100 10 1/1000))))))
1.832The simulated approach seems more intuitive to me with ~8% for a p-value of 5% and ~1.8% for a p-value of 1%.
P(bowl 1 | plain) = ( P(bowl 1) * P(plain | bowl 1) ) / P(plain)
= ( 0.5 * 30/40 ) / 50/80
= 0.6
(p/bayes 0.5 30/40 50/80)
0.6
P(94 bag | yellow) = ( P(94 bag) * P(yellow | 94 bag) ) / P(yellow)
P(94 bag) = 0.5
P(yellow) = 0.5 * 2 + 0.5 * 0.14
P(yellow | 94 bag) = 0.2
(p/bayes 0.5 0.2 (+ (* 0.5 0.2) (* 0.5 0.14)))
0.5882
P(twins) = 0.019
P(identical twins) = 0.002
P(fraternal twins) = 0.019 - 0.002 = 0.017
P(A|B) = P(A)*P(B|A)/P(B)
P(A|B)*P(B) = P(A)*P(B|A)
P(A|B)*P(B) = P(A and B)
P(A|B) = P(A and B)/P(B)
P(elvis was identical | twin brother) = P(elvis was identical) * P(twin brother | identical) / P(twin brother)
P(elvis was identical | twin brother) = P(elvis was identical and twin brother) / P(twin brother)
P(twin brother) = P(identical & twin brother) + P(fraternal & twin brother)
P(identical & twin brother) = 0.002 * 0.5
P(fraternal & twin brother) = 0.017 * 0.5 * 0.5
P(twin brother) = (0.002 * 0.5) + (0.017 * 0.5 * .05)
P(elvis was identical | twin brother) = ( 0.002 * 0.5 ) / ( (0.002 * 0.5) + (0.017 * 0.5 * 0.5) )
.19047
About 19%?